A. *^W"S of 



♦, --^ A m^<w A w,^ v ..... . 














.<°^ 



s^. *7^T» v*- ^ '*.To° .0" 














*°* 



v* * 











.0 4 



40. 




9 A^^ 






"* AT "*A • 
































* aV "^s •@Hl9> «? V 



« *J^fe!' *^K A* 



^ '^aws ** 





JP-V 










-W 






^ 



6°+ 






• ••' 






*••!•%% * <#sj^% V.^-V vv^rv 







G v V>. 












JO, • 








\5 



A' 

















^ 
*> 



.'dfeir. \/ ; m v^ 




^c> 'o . » * A <. ♦777T* .G v ^ 



,i°* 







^ -»• 







<^ a0 v •!••* > 



••JS&-X ftitifcS J*\-&\ * 4&k ° 




J°+ .\ 



*°"*v v 



^°- 







C27\ 



BUREAU OF MINES 
INFORMATION CIRCULAR/1990 



BOMSPS: Bureau of Mines Signal 
Processing Software-Concepts, 
Expressions, and Tutorial 



By Michael J. Friedel 




BUREAU OF MINES 
1910-1990 

THE MINERALS SOURCE 



Mission: As the Nation's principal conservation 
agency, the Department of the Interior has respon- 
sibility for most of our nationally-owned public 
lands and natural and cultural resources. This 
includes fostering wise use of our land and water 
resources, protecting our fish and wildlife, pre- 
serving the environmental and cultural values of 
our national parks and historical places, and pro- 
viding for the enjoyment of life through outdoor 
recreation. The Department assesses our energy 
and mineral resources and works to assure that 
their development is in the best interests of all 
our people. The Department also promotes the 
goals of the Take Pride in America campaign by 
encouraging stewardship and citizen responsibil- 
ity for the public lands and promoting citizen par- 
ticipation in their care. The Department also has 
a major responsibility for American Indian reser- 
vation communities and for people who live in 
Island Territories under U.S. Administration. 



{yS*^»- ^^£ V 



Information Circular/9242 



BOMSPS: Bureau of Mines Signal 
Processing Software-Concepts, 
Expressions, and Tutorial 



By Michael J. Friedel 



UNITED STATES DEPARTMENT OF THE INTERIOR 
Manuel Lujan, Jr., Secretary 

BUREAU OF MINES 
T S Ary, Director 



^ UP 



Library of Congress Cataloging in Publication Data: 



Friedel, Michael J. 

BOMSPS : Bureau of Mines signal processing software : concepts, expressions, 
and tutorial / by Michael J. Friedel. 

p. cm. - (Information circular / Bureau of Mines; 9242) 

Includes bibliographical references. 

Supt. of Docs, no.: I 28:27:9242. 

1. BOMSPS (computer program) 2. Seismic prospecting-Data processing. 
3. Signal processing-Digital techniques-Data processing. I. Title. II. Title: Bureau 
of Mines signal processing software. III. Series: Information circular (United States. 
Bureau of Mines); 9242. 

TN295.U4 [TN269] 622 s-dc20 [622M592'02855369] 89-600335 

CIP 



5 




CONTENTS 

Page 

Abstract 1 

Introduction 2 

Digital processing concepts and expressions 2 

Time series function 2 

Theoretical considerations 3 

Numerical considerations 5 

Digital filtering 9 

Single time series function 9 

Two time series functions 9 

Tutorial 14 

Program options 15 

Available time series functions 16 

Correlation modules 17 

(De) convolution module 17 

Filter module 18 

Fourier module 18 

Gain module 21 

Hilbert module 22 

Interpolate module 22 

Mute module 24 

Normalization module 25 

Plot module 26 

Polarity module 27 

Randomization module 27 

Scale module 28 

Shift module 28 

Smooth module 29 

Stack module 30 

Window module 30 

Write module 31 

Summary 32 

References 32 

Appendix A.-Sample data files 33 

Appendix B. -Mathematical relationship of Fourier transformation to Hilbert transformation 39 

Appendix C.-Symbols used in this report 40 

ILLUSTRATIONS 

1. Analog versus digital representation of time series function 3 

2. Flowchart depicting typical analog-to-digital conversion and data transfer schemes 4 

3. Complex representation of an arbitrary number and time series function 5 

4. Effect of sample interval on Nyquist frequency 6 

5. Aliasing phenomenon 7 

6. Frequency-domain representation of a boxcar function employing a fast Fourier transform 8 

7. Time-domain and frequency-domain representation of a Dirac delta function 9 

8. Schematic of three-component linear system 11 

9. Graphical representation of convolution property 11 

10. Graphical representation of cross-correlation property 13 

11. Application of linear systems analysis to band-pass filtering 13 

12. Frequency-domain filters 19 

13. Transient response to implementing frequency-domain filtering with regard to a Dirac delta function ... 20 

14. Screen view of trace attributes: original transient (input signal), quadrature, instantaneous amplitude 

(envelope), and instantaneous phase functions 23 

15. Interpolation of a digital time sequence 25 

16. Some time-domain operations as applied to a time series function 25 



ILLUSTRATIONS-Continued 



17. Screen view of Fourier attributes: transformation (real and imaginary), amplitude, and phase spectra 

18. Time-domain windowing cause, effect, and solution 



Page 

27 
29 



UNIT OF MEASURE ABBREVIATIONS USED IN THIS REPORT 


c/s 


cycle per second 


/is 


microsecond 


dB 


decibel 


rad 


radian 


Hz 


hertz 


rad/s 


radian per second 


kHz 


kilohertz 


s 


second 


MHz 


megahertz 


V 


volt 


ms 


millisecond 







BOMSPS: BUREAU OF MINES SIGNAL PROCESSING 
SOFTWARE-CONCEPTS, EXPRESSIONS, AND TUTORIAL 



By Michael J. Friedel 1 



ABSTRACT 

This information circular describes a comprehensive digital signal processing system, developed by 
the Bureau of Mines, called BOMSPS. Operations such as gain control, muting, normalization, polarity 
change, randomization, scaling, shifting, stacking, smoothing, and windowing account for the basic time- 
domain operations available. Higher order time-domain enhancement capabilities involve convolution, 
deconvolution, autocorrelation, and cross correlation. Frequency-domain attributes such as the real 
transform, imaginary transform, relative amplitude, and power and phase spectra can be determined. 
Frequency-domain filter capabilities permit the investigator to employ a half-sine or cosine window, 
notch band-pass, or high-cut or low-cut filter, while interactively defining the rolloff of either a linear 
or cosine taper. Calculation of related trace attributes such as quadrature (Hilbert transform), 
instantaneous amplitude, and instantaneous phase are also available. This program should be a welcome 
addition for those mining and environmental geophysicists, researchers, and engineers interested in 
processing signals. 



'Geophysicist, Twin Cities Research Center, Bureau of Mines, Minneapolis, MN. 



INTRODUCTION 



Both mining and petroleum geophysicists utilize seismic 
signals to characterize rock masses for delineating sub- 
surface stratigraphy and structural conditions. However, 
the vast amount of information available in project data 
(i.e., body, head, surface, airwave modes, and noise) often 
precludes any attempt at interpretation unless these data 
are enhanced. New, innovative processing technology 
increases the quality of these signals for analysis and in- 
terpretation. Seismic signal processing software, however, 
is written mostly for oil exploration companies. Conse- 
quently, commercially available programs (1) are often 
limited to a single application, (2) are typically restricted 
for use with specific hardware, e.g., a digital oscilloscope 
or a mainframe computer, (3) use computer codes that are 
usually proprietary and that restrict in-house modification 
to expand capabilities, and (4) are typically hard to in- 
terface, poorly documented, expensive, and support 
dependent. 

As part of its mission to assist the mining industry in 
the advancement of technology, the Bureau of Mines has 
developed BOMSPS (Bureau of Mines Signal Processing 
Software) to meet the need for a signal processing 
program. BOMSPS is a comprehensive software package 
that provides for improved signal quality to assist in the 
interpretation of seismic data collected in mining oper- 
ations. The software is capable of analyzing various time- 
and frequency-domain attributes and possesses filtering 
capabilities. Postdigital enhancement by BOMSPS to 
band-limit or reject selected information provides greater 
resolution for necessary interpretations. This software is 
an inexpensive tool for the geophysicist, engineer, and/or 



researcher performing seismic surveys in mining and en- 
vironmental operations. 

BOMSPS is a comprehensive digital processing inter- 
face system, written in FORTRAN 77 language for use on 
most personal computers. 2 Utilizing a fast microcomputer 
with relatively high volume on-line disk storage is an ef- 
ficient way to process the time series data. Software mod- 
ules incorporated into BOMSPS typically provide a single 
solution to a time- or frequency-domain operation. 

The primary objective of this report is to describe in 
detail the signal processing concepts and expressions in- 
corporated into BOMSPS. Although the text briefly dis- 
cusses the mathematical concept of the one-dimensional 
Fourier transformation, it omits proofs, which can be 
found in many textbooks. The review highlights associated 
Fourier properties and their use in developing mathemat- 
ical expressions used as digital operators. 

The second objective is to provide a user's document 
for the implementation of BOMSPS. Easy-to-follow 
menus and features that heretofore may have been found 
only in expensive code for mainframe computers, such as 
those codes available to the oil industry, are incorporated. 
Menu information and messages appearing on the screen 
during execution of BOMSPS are explained throughout the 
text. Sample data files and output are incorporated for 
ease of use by technical as well as nontechnical personnel. 
It is expected that the reader will gain a greater knowledge 
of both the advantages and limitations of these phases of 
digital processing, which are made available in BOMSPS 
for transient analysis. 



DIGITAL PROCESSING CONCEPTS AND EXPRESSIONS 



TIME SERIES FUNCTION 

A considerable amount of information is contained in 
a convoluted time series function (also called a transient, 
trace, signal, or waveform). The capacity of a time func- 
tion to carry useful information is related to a product of 
amplitude, length of the signal, and frequency bandwidth. 
Any continuous measurement of some attribute with re- 
spect to time is defined as an analog signal. Specifically, 
the independent variable is defined for all time, or s(t). 
Conversely, if the signal is defined by some number (N) of 
discrete points, it represents the analog signal's digital 
counterpart. In this sense, the digitized transient may 
represent a discrete version of the analog signal. 

In computer analysis of time functions, it is necessary 
to perform operations on the digital record. Figure 1 
depicts the analog and digital representations of an ar- 
bitrary time series function (SINUSOID.DAT, appendix 
A). The ordinate values are relative and normalized, while 
the abscissa values represent time units (seconds). The 



continuous line in figure 1A connotes a continuous 
sampling of the transient (e.g., by means of strip chart re- 
corder). In figure IB, the time function is defined by 
100 discrete points sampled at 1-ms intervals. Where 
discrete data points are depicted in subsequent figures, a 
line is drawn connecting the data points to visually high- 
light the signal under investigation. In order to illustrate 
other concepts, some figures omit sample points. The 
transients illustrated, however, are not to be assumed to 
be analog, since all numerical operations must be con- 
ducted on digital data. 



A source code listing is available upon request from M. J. Friedel, 
Twin Cities Research Center. The Bureau of Mines expressly declares 
that there are no warranties expressed or implied which apply to the 
software described herein. By acceptance and use of such software, 
which is conveyed to the user without consideration by the Bureau of 
Mines, the user hereof expressly waives any and all claims for damage 
and/or suits for or by reason of personal injury, or property damage, 
including special, consequential or other similar damages arising out of 
or in any way connected with the use of the software described herein. 



-1.0 




LU 
> 

I— 
< 

_J 
LU 



-.2 - 



.6 " 



10 



□ 


— i 1 1 i 




D 

n 


KEY 


B 


° m 

on o Sample point 




■°0 Q 


□ 




D 


D ^""*w 






B ^^. 




a a 


n jF X 




D D 


■ 1 \ A . 


" 


a 


□ a W_ Jo" 




D 


a D 






o o 






» o 




D a 


V 


- 


an 






■B 


1 1 1 1 





0.02 



0.04 



0.06 



0.08 



0.1 



TIME, s 

Figure 1. -Analog versus digital representation of time series 
function. A, Analog representation, continuous sampling; 
B, digitized representation, 100 discrete data points sampled at 
1-ms intervals. 



The analog-to-digital (A-D) conversion is handled in a 
number of ways. The simplest method is manually picking 
discrete amplitude values as a function of time. For mul- 
tiple signals exhibiting broadband spectrum and a relatively 
long time duration, discretization becomes a laborious 
task. While graphics tablets can be used, it is more com- 
mon to make A-D conversions with commercially available 
products such as converters, digital oscilloscopes, or seis- 
mographs. Also, A-D converters can be installed directly 
in a computer. As a result, analog data can be acquired 
and converted directly at the field or laboratory site. As 
appealing as this sounds, this type of conversion is not 
always possible. E.g., waveforms often are recorded at a 
field or laboratory site on either paper or magnetic tape. 
If the A-D process does not utilize an acquisition- 
processing computer, it is often necessary to digitize and 
store the data until they can be transferred to the central 
processing computer. This process necessitates additional 
software and hardware. Ultimately, there are additional 
steps required prior to the actual processing of sampled 



time-sequence data. A flowchart depicting typical 
acquisition-transfer schemes is given in figure 2. There 
are a number of vendors that can supply hardware and 
software for the digitization and transfer process. In the 
following discussions, it is assumed that the data are in 
ready-to-process digital format. 

THEORETICAL CONSIDERATIONS 

Given any physically realizable signal s(t) defined -°° to 
a>, there exists a Fourier transform pair representative of 
the signal in either time or frequency (I). 3 By definition, 
the one-dimensional transform pair includes a forward and 
an inverse expression. The Fourier transform pair pro- 
vides the fundamental vehicle by which BOMSPS maps 
information from one domain to the other (time to fre- 
quency, or frequency to time). The forward transform is 
defined by the following expression: 



S(w) = 



s(t) exp(jo?t)dt o s(t), 



(1) 



where S(u>) = frequency-domain equivalent of s(t), 

s(t) = time series function, 

u> = angular frequency (2nf), rad/s, 

f = frequency = 1/period, c/s or Hz, 

t = time, s, 

J - i-lf, 

n ~ 3.14159, rad, 

and <f> = 27rt = phase. 

The inverse transform is given by 



s(t) = 



S(w) exp(-ju>t)dw o S(w). 



(2) 



The Fourier transformation is designated by a two- 
headed arrow. 

In general, the Fourier transform is a complex quantity, 

S(«) = R(w) + jlm(u>) = | SO) | expGwt), (3) 

where R(^), I(w) = real and imaginary parts of the 
Fourier spectrum, respectively. 



Italic numbers in parentheses refer to items in the list of references 
preceding the appendixes. 




Transducer(s) 


































FM tape 








1 








Datagraph 














ir 










Seismograph 


Digital 
oscilloscope 








Graphics 
tablet 


















1 


r 






r \ 
Microcomputer 

v ) 


■^ 


r 











Central computing facility 



Analog (A) 
signal 

A-D 
conversion 

Digital (D) 
signal 



Figure 2.-Flowchart depicting typical analog-to-digital conversion and data transfer schemes. 



| S(u>) | represents the amplitude spectrum of the signal s(t) 
and is defined as 



|S(c*)| =[r 2 (u) + £(«)] 2 - (4) 



The phase spectrum (or angle) is defined as 



<f>(u>) = tan flm(w)/R(w)| . (5) 



These transform characteristics are depicted in figure 3. 
In figure 3A, a two-dimensional representation of a time 
point is shown. The corresponding phase angle (initial 
angle) and vector magnitude (length of vector) are given. 
The amplitude is a resultant vector between the real and 
imaginary axes. The angle that this vector makes with the 
real axis is defined as the phase. In figure 3B, this idea is 
extended to include the transformation of a continuous 
time function. The relationships that govern these con- 
cepts are 



and 



R(w) = |S(w)|cos(wt) 
Im(w) = -j|S(w)|sin(wt). 



(6) 
(7) 



3 

E 



LU 

z: 
o 

y: 
o 
u 

>- 

CL 
< 

CD 
< 




Rotating 
vectors 



w (rad/s) 




REAL COMPONENT [R(w)] 



Figure 3.-Complex representation of an arbitrary number and 
time series function. 



As a consequence, the expanded form of the Fourier trans- 
form can be written as follows: 



S(w) = | S(w) | cos(wt)dt - j | S(w) | sin(wt)dt. (8) 

Since cos(wt) = cos(-wt), it immediately follows that the 
real part of the transformation is R(-w) = R(w). There- 
fore, the real component is an even function. Similarly, 
sin(-a)t) = -sin(o)t), which implies that I m (w) = -Im(w). 
Hence, one can assume that the imaginary components of 
the transformation represent an odd function. Since an 
odd function is equal to zero between symmetric limits 
(Fourier integral), the amplitude function presents itself as 
an even function. A similar analysis establishes that the 
phase spectrum is an odd function. 



NUMERICAL CONSIDERATIONS 

The wide range of problems susceptible to analysis by 
the Fourier transform (analytic solution) led to the de- 
velopment of the discrete Fourier transform (DFT) (2-3). 
The DFT is the numerical equivalent of the analytic ex- 
pression (equation 8) shown above. The DFT formulation 
is based on the relationship 



S(f) = N sVt i )exp(j27rf k t i )(t i + 1 -t i ), 
k 1 = 



(9) 



where 



S(f) = frequency-domain representation of an- 
alog signal. 



This relation implies that there are N discrete data points 
sampled at an even spacing, hence, the name "discrete 
Fourier transform." In contrast to the analytic expression, 
the DFT operates on a signal with finite limits. Thus, it 
must be assumed that outside the given time window 
periodicity exists. This does not imply that the signal s(t) 
is an even function. 

Adequate sampling rates are necessary to preserve the 
integrity of the frequency content of a signal. Based on 
the Shannon sampling theorem (i), the time increment 
(At) between discrete points is defined as equal to the 
inverse of two times the maximum frequency (f ) within 
the wave train. 



At = l/2f Q . 



(10) 



This is an important consideration in determining the 
Nyquist frequency (f ), also known as the folding or 
aliasing frequency, defined as 



f nq = V2At. 



(11) 



It is the Nyquist frequency that is the maximum that can 
be examined within the wave train. A common but erro- 
neous belief is that additional sampling (increased sample 
rate) will always provide greater resolution. This is true 
only when detailed information (higher relative frequen- 
cies) is actually present in the signal. 

Once the sample interval (or sampling rate) is suffi- 
cient to define the maximum spectral component in the 
waveform, no additional information can be obtained. 
Figure 4 depicts the effect of perturbing sample interval 
on the Nyquist frequency. The signal under investigation 
represents an arbitrary sinusoid (SINUSOID.DAT, Ap- 
pendix A). In figure 4A, the sample interval is 1 ms 
(1 kHz), resulting in a Nyquist frequency of 500 Hz. By 
decreasing the sample interval to 125-/zs (8.2 MHz), the 
Nyquist frequency is effectively shifted to 4,000 Hz. 



1.0 
.8 

.6 



LU 

Q 


.2 


Z) 




1— 




_l 




Ql 




7" 




< 





LU 




> 




1— 
< 


1.0 


_l 




LU 




CH 


.8 




.6 




100 



200 



300 



400 



500 



1 1 1 1 


i 


Time 3erie3 function: 


SINUS0ID.DAT, appendix A 


Sample interval 


= 125U3 


1 Nyqui3t frequency 


J = 4,000 Hz 


Maximum frequency 


1 = 500 Hz 



800 1,600 2,400 3,200 4,000 

FREQUENCY, Hz 

Figure 4.-Effect of sample interval on Nyquist frequency on 
a time series function (SINUSOID.DAT, appendix A). 

Decreasing the sample interval (increasing the sample 
rate) beyond this will extend the Nyquist frequency into a 
region devoid of signal energy. In addition, excessive sam- 
pling may contribute to a degradation in the signal-to-noise 
(S-N) ratio and may require an additional amount of array 
storage. In figure 4B, there is no additional energy above 
about 500 Hz. 

If the sampled sequence contains frequencies higher 
than one-half of the Nyquist frequency, those higher fre- 
quencies will be lost. The sample frequency (Af) is defined 
by 

Af = 2f nq /NP2, (12) 

where NP2 = total number of discrete points associated 
with time sequence. 



With aliasing of a signal, high frequencies of the original 
spectrum are "folded" back on top of valid lower 
frequencies. Hence, there is (1) a loss of high-frequency 
energy and (2) distortion of the low-frequency energy. 
While aliasing prohibits a meaningful frequency-domain 
interpretation, it completely destroys the originally sampled 
time function. Figure 5 shows the phenomenon of aliasing. 
Figure 5A depicts the analog time series function and 
corresponding frequency spectrum. The erroneously sam- 
pled signal and corresponding frequency spectrum is shown 
in figure 5B. Folding of the aliased frequencies is shown 
in 5C. Inspection of the aliased signal reveals little sim- 
ilarity to its original counterpart. The information lost 
during poor sampling cannot be recovered. 

A common problem that occurs during the data acqui- 
sition phase is in selecting an aliasing (high-cut) filter that 
is approximately equal to the Nyquist frequency. A given 
high-cut filter does not eradicate all frequencies above its 
designated value. In practice, the alias filters have cutoff 
slopes, beyond the designated frequency, of roughly 70 dB 
per octave. The reduction in frequency is about 1/2000 
in a doubling of frequency, e.g., from 125 to 250 Hz for a 
2-ms sampling interval. Clearly, proper implementation of 
an antialiasing filter would require familiarity with its 
response character. The aliasing filter used should be one 
that has no appreciable response at the Nyquist frequency. 
Other practical considerations involving knowledge of the 
Nyquist frequency include applications of frequency- 
domain filtering and linear systems analysis. Both of these 
topics are addressed later in this report, in the section 
"Two Time Series Functions," under "Digital Filtering." 

To determine the amplitude of N separate sinusoids, 
the computational time is found to be proportional to the 
number of multiplications (N). Since the computational 
time associated with the discrete transformation can be 
prohibitive for most applications (in terms of time and 
expense), the fast Fourier transform (FFT) was devised. 
The FFT is a numerical algorithm that reduces the com- 
putational speed to a time proportional to N log N. For 
N = 2, the FFT algorithm factors an N x N matrix, so that 
each of the factored matrices minimizes the number of 
multiplications and additions. 

There are many FFT algorithms available on the soft- 
ware market. BOMSPS employs a three-dimensional FFT 
developed by Robinson (4), but others could be used. 
Only one dimension, however, is actively required by 
BOMSPS. Future signal processing capabilities may be 
extended to either two or three dimensions. The FFT 
requires that the number of discrete points to be operated 
on is a power of 2. Most external data files contain some 
arbitrary number of data points. In order to comply with 
the power-of-2 number requirement, external data files are 
automatically padded with additional data points (of zero 
amplitude) until the power-of-2 requirement is accom- 
modated. The padding of data files is handled by the 
NPTWO module during execution of BOMSPS. 




Z) 



Q_ 



< 

> 



< 



CL 



-VWVyi 



<^> 



Nyqui3t 
frequency 



B 



££\ 



Sample point 




TIME 



FREQUENCY 



Figure 5,-Aliasing phenomenon. A, Analog time series function and corresponding amplitude spectrum; B, proposed time sampling 
and corresponding amplitude spectrum; C, resulting aliased signal and corresponding amplitude spectrum. The shading in C indicates 
the "folding" of higher frequencies back on valid frequencies. The two-headed arrow represents the transformation process. 



The probability that the originally sampled time se- 
quence is an odd function is great, particularly if the data 
are not a contrived function. The FFT performs compu- 
tations strictly in positive space. Specifically, all negative 
frequency values are mapped as positive values. Corre- 
sponding arrays are characteristic of increasing positive 
frequency values from Hz to the Nyquist frequency 
(number of N points divided by 2, NP2 points). Beyond 
this point, the maximum negative frequencies begin at the 
Nyquist frequency plus 1 frequency increment (NP2 + 1 
points) to 2 times the Nyquist frequency. Hence, there 
exits symmetry about the Nyquist frequency and not about 
Hz. 



The characteristic amplitude spectrum (Sine Function) 
of a boxcar function with a period of 40 ms, sampled ev- 
ery 1 ms, and of unit amplitude (BOXCAR.DAT, appen- 
dix A), is displayed in figure 6. This spectrum is indica- 
tive of the manner in which FFT operates (in positive 
space). Relative amplitude values stored in the compu- 
tational array appear to begin at Hz and extend out to 
1,000 Hz (1 kHz). In reality, those frequency values that 
are depicted beyond the folding frequency are negative 
values (as previously described). The true spectral equiv- 
alences are depicted along the lower edge of the figure. 
An understanding of FFT operation and storage of calcu- 
lated values is important in considering an extension 



200 



400 



600 



800 1,000 



Q 

Z) 



YL 
< 

LU 
> 

I— 
< 

_l 
LU 



1.0 



.4 " 



.2 " 




200 



400 



-400 



FREQUENCY. Hz 



Figure 6.-Frequency-domain representation (Sine function) of 
a boxcar function (BOXCAR.DAT, appendix A) employing a fast 
Fourier transform that operates in positive space. The top scale 
indicates computer storage of the spectrum, while the bottom 
portrays the actual spectrum. 



Clearly, the Hilbert transform must be numerically deter- 
mined prior to the calculation of the envelope. BOMSPS 
makes use of a Fourier-Hilbert transform relationship, in 
order to calculate the quadrature function. The mathe- 
matical relationship utilized for calculating the quadrature 
function (and indirectly the remaining Hilbert attributes) 
is outlined in appendix B. 

Interpretive use of the instantaneous amplitude is based 
on its relatively insensitive nature with respect to time. 
Specifically, the instantaneous amplitude is less sensitive 
to local variations caused by noise than the actual corre- 
sponding transient under investigation. Since relatively 
narrow band signals tend to oscillate and decay slowly, 
the apparent time of the transient response (and corre- 
sponding envelope) is extended. 

A display of the phase angle as a function of time is 
defined as the instantaneous phase. The instantaneous 
phase physically amounts to removing all amplitude infor- 
mation. BOMSPS utilizes the following mathematical 
relationship to calculate the instantaneous phase &(t): 



^i(t) = tan 



(Ws(t)]" 1 . 



(14) 



of the existing FORTRAN code. Currently, those spectral 
functions written to external data files are characteristic 
of only one-half of any given spectrum (0 Hz to the 
Nyquist frequency). If a need arises to write both positive 
and negative spectra, the user can modify the WRITE 
statements (implied DO loops) to go from 1 to NP2, 
instead of looping from 1 to NP2/2. 

Time series attributes, such as the instantaneous am- 
plitude (envelope), instantaneous phase, instantaneous 
frequency, and derivative, are calculated by making use of 
the quadrature function, or Hilbert transform of the orig- 
inal signal. Hence, these functions fall under the heading 
of "Hilbert attributes." E.g., the expression related to the 
envelope e(t) (instantaneous amplitude) of a time function 
is 



e(t) = [s(t) + s\t)j , (13) 



The instantaneous phase angle calculated above will wrap 
around the polar axis of the sampled time series. Between 
any two successive discrete points, the amount of phase 
shift can be determined. If the phase shift between two 
sample points is multiplied by the sample rate and divided 
by 360 (degrees), one will arrive at the instantaneous 
frequency f 4 (t). While this is currently not incorporated 
into the HILBERT module, the extension to include this 



is 



fj(t) = d^(t)/dt. 



(15) 



Similarly, the derivative is not currently an option avail- 
able to the user. However, the mathematical expression 
given below can be implemented numerically, based on 
previously defined functions. 



ds(t)/dt = -jwS(w). 



(16) 



where s(t) = the input time function, 



and s(t) = Hilbert transform of this signal. 



DIGITAL FILTERING 



SINGLE TIME SERIES FUNCTION 

Digital filtering provides methods to discriminate 
against certain transient attributes for data enhancement. 
The filter (a function of time or frequency) is numerically 
multiplied by the time- or frequency-domain representation 
of the signal under investigation. The simplest digital filter 
involves multiplication of a given time transient by a con- 
stant (i.e., scaling). The application of both frequency- 
and time-domain filters makes use of the Fourier trans- 
form scaling property. The following dependence exists 
for scaled time- or frequency-domain functions. 



b(ct) ** 1/c B(w/c). 
By symmetry relations, 



1/c b(t/c) ~ B(cw), 



(17) 



(18) 



where b(t), B(w) 



time- and frequency-domain equiv- 
alents for the boxcar function, 
respectively, 



and 



c = arbitrary constant. 



By inspection it can be determined that time-scale ex- 
pansion corresponds to frequency-domain compression. 
Similarly, as the frequency scale expands, the corre- 
sponding time scale is compressed. This might have been 
expected based on the symmetry property of Fourier trans- 
forms. As the time or frequency scale expands, the am- 
plitude of corresponding transient quantity (frequency or 
time) increases, to keep the area constant. The concept 
of Fourier scaling can best be demonstrated by employing 
the Dirac delta function during this process. Figure 7 
depicts the transformation of the Dirac delta function. 
Note that there are no side lobes present (as opposed to 
the Sine function, figure 6). More specifically, the Fourier 
transform of the Dirac delta function contains the com- 
plete spectrum, all of equal amplitude. The relative am- 
plitude is equivalent to that of its time-domain represen- 
tation. This property of the Dirac delta function is used in 
the discussion of convolution in the following section. 

TWO TIME SERIES FUNCTIONS 

Up to this point, all filtering concepts are based on the 
scaling property of the Fourier transform. If the functions 
a(t) and b(t) have Fourier transforms A(w) and B(w), 
respectively, then the sum a(t) + b(t) has the transform 
equivalent of A(w) + B(w). The transform pair can be 
written as follows: 



This Fourier property establishes the applicability of ana- 
lyzing a linear system by Fourier transformation. Using 
Fourier tranformation, the linear systems approach can be 
employed as a supplemental approach to digital filtering. 
Concepts based on this approach include cross correlation, 
autocorrelation, convolution, and deconvolution. The 
investigator may choose to invoke any one of these numer- 
ical operations for filtering applications, in either the time 
or frequency domain. Upon successful implementation 
of any one of these operations, the user may incorporate 
any (or all) of the remaining menu options, to facilitate 



2 

1.6 

> 

Q I 2 



< 



.8 

.4 



I ' 1 

■ ■ 


— r ■■ ■ i r I 

<l 



0.08 0.16 0.24 0.32 0.4 

TIME.s 




a(t) + b(t) ** A(w) + B(w). 



(19) 






Cl. 

H 

< 

LU 
> 

I— 
< 

_l 
LU 

CL 



FREQUENCY, Hz 

Figure 7. -Time-domain (A) and frequency-domain 
(B) representation of a Dirac delta function shifted 200 ms 
(DIRAC.DAT, appendix). The spectrum is darkened in to illustrate 
the transformation. 



10 



additional data enhancement. The following mathematical 
relations are utilized in the computational algorithms em- 
ployed for convolution and/or deconvolution of a linear 
system: 

Convolution: 

f(t), i(t), o(t) - F(co), I(u), O(u), (20) 

f(t) * i(t) = o(t) - F(u) 1(c) = O(o>), (21) 

f 



f(t) * i(t) = 



(r?)i(t-r?)df?; 



(22) 



Deconvolution: 



where 



f(t) - F(o>) = 0(a>) / I(a>), (23) 

i(t) ~ I(w) = O(w) / F(w), (24) 

* = time-domain convolution operator, 
f(t) = forcing function, 
i(t) = impulse response function, 
o(t) = transient response function, 



and 



F(u>), I(w), O(w) = transform equivalents of forcing, 
impulse, and transient response 
functions, respectively. 

The frequency-domain equivalent of the impulse response 
function is often called the transfer function; however, the 
other functions are simply denoted as the forcing and 
response function transform equivalents. A schematic of 
this three-component linear system is shown in figure 8. 
Since autocorrelation is a special case of cross correlation, 
and the operations of convolution and deconvolution are 
related, the following text is broken into two discussions. 
In general, the time-domain convolution process 
(TDCON module) results from a continuous multiple-add 
sequence. The responses for each of a series of impulses 
are summed at the appropriate spacing defined by the 
sample interval. Numerically speaking, the signal to be 
shifted (operator) is first folded, or reversed end for end. 
Next, the folded signal is placed to the left of the digitized 
input transient. The operator is shifted right one position. 
Each of the coincident values of the operator and signal 
are then multiplied. All the products are then summed to 



yield one output sample corresponding to that time value. 
The operator is then shifted right one position and the 
operation repeated (until all positions have been occu- 
pied), resulting in the response function. The convolution 
operation is depicted in figure 9. 

The frequency-domain convolution (FDCONV module) 
approach requires that both the signal and operator be 
transformed using the FFT. Each coincident value is then 
multiplied. The function resulting from the product of 
these two functions is then inverted back in time, in order 
to yield the fdtered transient response. Since multiplica- 
tion in the transform domain is equivalent to the multiple- 
add sequence conducted in the time domain, both transient 
response functions are theoretically identical. However, 
numerical algorithms can reproduce poor results depen- 
dent on certain characteristics of the sampled time se- 
quence under investigation. In a similar manner, the 
deconvolution process (FDDCON module) requires that 
both functions be transformed. In this case, division of 
coincident values occurs. Inverting the function established 
during the frequency-domain deconvolution results in 
producing either a forcing or an impulse function. 

The two mathematical operations convolution and cross 
correlation, account for the majority of basic signal pro- 
cessing. As it did with convolution and deconvolution, 
BOMSPS provides the user with alternatives for imple- 
menting cross correlation in either time or frequency. The 
implementation of cross correlation suggests that there 
exists a linear relationship between the forcing, impulse 
response, and transient response functions. In general, the 
correlation function represents a graph of the similarity 
between two waveforms. If a relative time shift(s) is al- 
lowed between two waveforms, then a correlation function 
can be described. E.g., the correlation of a pure sinusoid 
(monotonic frequency) will result in extracting that same 
frequency domain, while rejecting the others. Clearly, such 
a property is advantageous to band-pass filtering. Mathe- 
matical relationships defining the cross correlation and 
related frequency-domain attributes of a signal are given 
by 



f(t)i(r + t)dt 



M> = ; -; * *(«) = F («)!(«), ( 2 5) 

c ( ( ( ^' 2 r 



f*(t)dt 



i Z (t)dt 



*r,(t) * * fi («) = l* fi (w)| exp (ja>*), (26) 



11 



fCO 



i(0 



o(0 



Time domain 




Impulse response 



Transient 
response 



F (oo) x I (w) = (w) Frequency domain 

Figure 8.-Schematic of three-component linear system (* denotes convolution, and * denotes multiplication). 



2 


1 






□ 





1 


2 


3 


2 


1 


1 


Sigr 




i 


r 






1 

4 
7 
8 
5 
2 
2 


4 

4 

4 

4 

4 

4 

4 

4 



















1 +0 














2+2 +0 














4 + 3 + 














6+2 +0 














4+ 1 +0 














2+ 1 +0 














2 + 0+0 



► Operator (filter, i(t)) 

Signal (forcing function, f(t)) 



Transient response (o(t)) 

Figure 9. -Graphical representation of convolution property. 



where ^n(t) 

f(0, i(t) 



correlation function, 

forcing and impulse response 
functions, respectively, and 



I(w) = frequency-domain impulse response 
function, 



F (to) - conjugate transformation of forcing $fj( w ) and jwtf 

function, 



cross-power sprectrum and cross- 
phase spectrum. 



12 



In performing cross correlation in the frequency-domain, 
the two time functions must first be transformed. Then 
the conjugate function of the signal to be correlated is 
determined based on the relationship given above. The 
correlation operator is determined by taking the inverse 
Fourier transformation of this product. The absolute value 
of the correlation operator represents the cross-power 
spectrum. These functions are synonymous with the 
Fourier attributes (amplitude and phase spectra) defined 
for a single time function. However, these Fourier at- 
tributes now include the prefix "cross," which is indicative 
of a correlation operation. 

It is imperative that the correlation function be specified 
at the onset of the calculation, which is not true for the 
convolution process. To ensure that the proper operator 
is specified for cross-correlation calculations, BOMSPS 
requires linear functions to be defined when any linear 
operation is performed. The only operational time-domain 
difference between convolution and cross correlation is 
that the cross-correlation operator is not folded. Coin- 
cident values are multiplied and summed to yield a single 
data point. The operator is then moved one position and 
the process repeated. This process is displayed in fig- 
ure 10. In this case, the forcing function (input) and op- 
erator are characterized by the same functions utilized in 
the convolution example. Specifically, the forcing function 
used is the sinusoid (SINUSOID.DAT, appendix A), while 
the impulse response function is an arbitrary function. 

In general, the correlation process serves as a band-pass 
filter, passing common frequency components between the 
signal (forcing function) and filter (operator). All other 
frequencies are effectively rejected. By cross-correlating 
the received signal (plus noise) against a known input 
signal, a cross-correlated peak is obtained at the lag time 
corresponding to the arrival time of the signal. This fact 
is important in determining the apparent velocity between 
these two points. While correlation remains insensitive to 
amplitude differences between waveforms, it is sensitive to 
stretching and/or compression of a waveform. 

The result of employing a Dirac delta function as the 
filter (impulse response-transfer function) in the previ- 
ously mentioned processing schemes is depicted in fig- 
ure 11. From the diagram, it is seen that the transient 
response function is identical to the forcing function 
(SIGNAL.DAT). This phenomenon will occur regardless 
of the operation utilized (convolution or correlation). This 
is true, however, only when a Dirac delta function centered 
at zero time is employed. While an impulse response 
characterized by a shifted Dirac delta function will result 
in the same transient response shape, time positioning will 
be different and dependent on the operation employed. If 
the operator used is not a Dirac delta function, the re- 
sulting transient response functions will take on a totally 
different character, each dependent on the operation em- 
ployed. One might have deduced this based on the scaling 
property of the Fourier transform (time-domain compres- 
sion represents complete frequency-domain expansion). 
Note that the transformed Dirac delta function (transfer 



function) results in a spectrum of constant unit amplitude 
of infinite bandwidth. In this case, the transfer function 
becomes an idealized band-pass filter, known as an all-pass 
filter. 

Another band-pass filter, useful for detecting hidden 
periodicities in a noisy waveform, is autocorrelation. The 
related frequency-domain information is the power density 
function. These functions are computed based on the 
following expression: 



^KO^fK") = l p (")l 



(27) 



where ^ ff (t) = autocorrelation function, 

and ^ff(w) = its frequency-domain representation. 

Autocorrelation is a special case of cross correlation. For 
this reason, the time- and frequency-domain processes of 
autocorrelation are identical to those discussed for cross 
correlation. However, now the forcing function is cor- 
related against itself (operator). The operation of auto- 
correlation, in the frequency domain, results in squaring 
the amplitude spectrum and subtracting phase values from 
themselves. The process of discarding all of the phase 
information consolidates any amplitude information. The 
periodic components of frequency in the operator and in 
the signal under investigation reappear unchanged in the 
autocorrelation function. Their amplitudes, however, are 
altered, giving greater prominence to the larger compo- 
nents. Based on these facts, the resulting time function 
would be expected to be of zero phase and to be sym- 
metric about the origin. The response results in a graph 
of the similarity between the signal and a time-shifted 
version of itself (correlation operator) with its corre- 
sponding maximum amplitude at Hz. 

Since the process of autocorrelation results in dis- 
carding the phase information, one could not expect to 
recover the original time function. Autocorrelation, how- 
ever, is useful in revealing characteristics about phenomena 
associated with the recording of the time function. The 
fact that an autocorrelation function cannot distinguish 
between the signal and background noise having similar 
spectra reveals its utility. In the event that a time function 
is truly random, the autocorrelation function will display a 
zero-amplitude time function. Hence, the autocorrelation 
of a signal inundated with random, uncorrectable transient 
activity will result in a spectrum free from noise. Quali- 
tatively, one may interpret the autocorrelated function in 
terms of the relative bandwidth. E.g., narrower band 
signals will oscillate over longer periods with the amplitude 
decaying at a slower rate. Furthermore, based on the 
Fourier scaling property, a broader band signal would be 
expected to exhibit a narrower autocorrelation function. 
Inspection of the power spectrum can then yield quantita- 
tive information as to the actual bandwith of the signal. 
Also, the ratio of the autocorrelation of a signal to the 
autocorrelation of noise yields the signal-to-noise (S-N) 
ratio. 



13 






1 


2 





X 


X 


X 












1 


2 


3 


2 


1 


1 


Sig 






r ^ r 1 






8 
7 
6 
3 
1 




4 — 
4 — 

4 

i 

4 

4 


+ 2+6 














+ 3+4 














+ 2+4 














+ 1+2 














+ 1+0 














+ 0+0 

» * 



Operator (filter, i(t)) 
Signal (forcing function, f(t)) 



Transient response (o(t)) 

Figure 10. -Graphical representation of cross-correlation property. 



f(t) 



i(t) 



o(0 



Time domain 




F(w) x i (co) = o(w) Frequency domain 

Figure 1 1 .-Application of linear systems analysis to band-pass filtering (* denotes convolution, and * denotes multiplication). 



14 

TUTORIAL 

To execute the signal processing program, the following Upon successful execution, the header listed below appears 
command sequence is entered after the prompt sign ( > ) on the screen. 
appears on the screen. 



>BOMSPS 



* BUREAU OF MINES * 
** SIGNAL PROCESSING INTERACTIVE SOFTWARE ** 
*** (BOMSPS, 1988) *** 

*** developed by *** 

** BUREAU of MINES ** 

* TWIN CITIES RESEARCH CENTER 
GEOTECHNOLOGY personnel 



* GEOTECHNOLOGY personnel * 



The first user-interactive message seen on the screen is options available to the user. The following example is 
directed toward input-output (I-O). Scrolling through from the program, 
these I-O messages results in listing I-O format and/or 

DO YOU WISH TO WRITE REMARKS CONCERNING 
INPUT/OUTPUT FORMAT OR PROGRAM OPTIONS ? 

ENTER: 1.0 (= YES) 
2.0 (= NO) 
>1 

EXTERNAL DATA FILES 

INPUT: Time per point (seconds), Total # data points 
Data points (time.amplitude) 

F10.8,I5 

F15.5,lx,F15.5 

Note: Total # data points <2048 

OUTPUT: To WRITE functional values to any/all of 

the data files to scratch files, the OPERATOR 
values (in external file WRT.IN) must be set 
equal to "1.0." Conversely, a value of "2.0" 
will prohibit the transfer of data. 

OPERATOR SCRATCH FILE FUNCTION 

(1) FISGNL.DAT FILTERED time fct (real component) 

(2) IMSGNL.DAT FILTERED time fct (imag.component) 

(3) QUAD.DAT QUADRATURE fct 

(4) ENV.DAT ENVELOPE (inst. ampltude) 

(5) PHI.DAT INSTANTANEOUS phase 

(6) TRANSR.DAT TRANSFORM spectrum (real component) 

(7) TRANSI.DAT TRANSFORM spectrum (imag.component) 

(8) AMP.DAT AMPLITUDE spectrum 

(9) PDS.DAT POWER spectrum 



15 

(10) PHAS.DAT PHASE spectrum 

(11) DAMP.DAT DYNAMIC spectrum (amplitude) 

(12) DPDS.DAT DYNAMIC spectrum (power) 

OUTPUT: Functional data values are written out to 
scratch files in the following format: 
f 15.5, lx,f 15.5 In this manner, a paper 
copy of the data (or screen visual) can be 
made employing other commercial graphics 
programs 

EXIT: To EXIT program at any point TYPE:CNTL C 

Based on the information presented above, the first line of Depending on the machine capabilities, this may be 

any input data file must include the sample interval in increased or decreased. Sample input data files are given 

seconds (F10.8) and the total number of data pairs (15). in appendix A. 

Information following this line represents the 

time-amplitude pairs (F15.5,lx,F15.5) with time increasing PROGRAM OPTIONS 

in succession. The maximum number of data points 

capable of being analyzed is based on the array space. After the request to review the processing options is 

Currently, the time-equivalent arrays are set to 2048. acknowledged, the following outline appears on the screen: 

AVAILABLE OPTIONS: 

A. SIGNAL PROCESSING 

1. SINGLE TIME SERIES FUNCTION 

a. SELECT from one of five standard signals 
or read in an external signal 

b. INTERPOLATE a random (or even) digitized 
data set to "new" even time interval 

c. SCALE time or amplitude by some constant 

d. STACK arbitrary number of waveforms 

e. GAIN control, convert arbitrary to true voltage 
values (and/or voltage values to velocities) 

f. MUTE any part of the input signal 

g. SHAPE input signal with either a Hanning or 
cosine weighting function, 3-pt smoothing 

h. RANDOMIZE (add random noise to input signal) 
i. TIME-SHIFT, translate function along time axis 
j. PHASE, modify signal phase (any amnt in radians) 
k. FILTER signal employing 60-Hz notch, bandpass, 

low or high cut with cosine or linear tapers 
1. CALCULATE any/ all of the following functions 



signal (real part) 
signal (imaginary part) 
cumulative time-domain energy 
envelope (instantaneous amplitude) 
quadrature (Hilbert transform of signal) 
instantaneous phase 
instantaneous frequency 

* derivative 

* transform spectrum 

* amplitude spectrum 

* phase spectrum 

* autocorrelation 

* power density spectrum 



16 



2. TWO TIME SERIES FUNCTIONS 

a. READ input signals from external files 

b. CALCULATE any/all of the following functions 
** cross correlation (time or freq. domain) 

* cross-power spectrum 

* cross-phase spectrum 

** (de) convolution (time or freq. domain) 

* transform spectrum 

* amplitude spectrum 

* phase spectrum 
** envelope 

** quadrature function 

** instantaneous phase 

** instantaneous frequency 

** derivative 



AVAILABLE TIME SERIES FUNCTIONS 

Since BOMSPS is written to be utilized on an instruc- 
tional, as well as a scientific, level, the user has access to 
various external time functions. The first menu to be 
displayed on the screen lists six available time functions for 
data analysis: 

MENU: AVAILABLE TIME SERIES FUNCTIONS 

1. BOXCAR (boxcar.dat) 

2. BOXCAR, shifted 100 ms . . (boxcarl.dat) 

3. DIRAC (dirac.dat) 

4. DIRAC, shifted 100 ms . . . (diracl.dat) 

5. EXTERNAL (filename.dat) 

6. SINUSOIDAL (sinusoid.dat) 

ENTER: time series function (filename.dat) = ? 

In addition to having the possibility of reading external 
data files, the user can investigate the response to a box- 
car, shifted boxcar, Dirac, shifted Dirac, frequency- 
modulated, and/or sinusoidal function. Each of these 
functions is characteristically sampled at 1-ms intervals. 
With the exception of externally created data files, the 
remaining functions exist in permanent data files (file 
names are given in parentheses above). The advantage of 
any numerical model is that it can be modified by simply 
changing a few values (i.e., sample interval, amplitude, 
etc.) Hence, the user may edit the file-perturbing values 
for extended analysis. In order for BOMSPS to read the 
appropriate data file, the file name is entered complete 
with the extension (filename.dat). In the example below, 
the data file containing an arbitrary sinusoid is selected by 
entering 

> SINUSOID.DAT 



For some computer systems the extension (.DAT) is in- 
cluded by default and does not need to be typed in. If the 
user wishes to investigate a field record (external option), 
the complete data file name must be typed. After the 
program has read a data file, the main processing menu 
appears on the screen as depicted below. 

*************************************************** 
SIGNAL PROCESSING MENU 



AUTOCORRELATE (time-domain). . 

CROSS CORRELATION 

(DE)CONVOLUTION 

DISPLAY 

EXIT PROGRAM 

FILTER (temporal) 

FOURIER-DOMAIN (attributes) 

GAIN/VELOCITY 

HILBERT-DOMAIN (attributes) 

INTERPOLATE 

MUTE 

NORMALIZE 

POLARITY 

RANDOMIZE 

RESET 

SCALE 

SHIFT 

STACK 

SMOOTHING 

WINDOW 

WRITE (values to screen) 

WRITE (values to scratch files) 



ENTER: 1.0 

2.0 

3.0 

4.0 

5.0 

6.0 

7.0 

8.0 

9.0 
10.0 
11.0 
12.0 
13.0 
14.0 
15.0 
16.0 
17.0 
18.0 
19.0 
20.0 
21.0 

22.0 

*************************************************** 



There are 22 operations available to the user, listed in 
alphabetical order in the menu. This order does not sug- 
gest any relative importance or endorse a particular pro- 
cessing scheme. 



17 



CORRELATION MODULES 

If cross correlation between a signal and operator 
(filter) is necessary, the investigator invokes option 2 
from the main menu. The following message appears, 
prompting the user for the file name in which the operator 
is located. 

ENTER: Cross correlation operator (filename.dat) = ? 
> dirac.dat 



(DE)CONVOLUTION MODULE 

The convolution-deconvolution operation is invoked by 
selecting option 3 from the signal processing menu. Im- 
mediately following this selection, the user is required to 
define the file name in which the operator (second time 
function, see below) is stored. 

ENTER: Time function # 2 (fdename.dat) = ? 
> diracl.dat 



Cross-correlating a signal with the unit spike results in the 
same output function as given in figure 11. This is the 
only condition in which these two operators will yield the 
same results. To determine the cross-power and cross- 
phase spectra, the user simply implements the FOURIER 
module (option 7). Note, however, that the screen display 
will retain the Fourier display titles associated with the 
single time transient. A special case of cross correlation 
is autocorrelation. The autocorrelation, or cross corre- 
lation of the signal with itself, can be conducted by using 
option 1. In this case, the signal is automatically assumed 
to be both the signal and operator. In a similar manner, 
the power and phase spectra can be determined by using 
the FOURIER module (option 7). 



The appropriate fde name is entered by typing the com- 
plete fde name and extension. In the example given, the 
file name used is DIRAC.DAT. This fde retains the Dirac 
delta function. BOMSPS routinely executes a comparative 
check of both time function sample intervals. If these 
values are not equivalent, the program will generate an 
error message, display the sample intervals for each trace, 
and then terminate. On the other hand, if the signals were 
sampled at concurrent rates, the user is requested to de- 
fine the linear system as defined below. 



PLEASE IDENTIFY THE INPUT TIME SERIES SIGNALS AS ONE 
OF THE CORRESPONDING FUNCTIONS IN THE LINEAR SYSTEM 



ENTER: 1.0 (= FORCING FUNCTION 

2.0 (= SYSTEM IMPULSE FUNCTION (OPERATOR) 
3.0 (= TRANSIENT RESPONSE FUNCTION 



As outlined above, two of the three linear functions must 
be defined (in time) as either forcing, system impulse 
(operator), or transient response (output) functions. In 
order to define the linear system, the user must identify 
each of the time functions as a component in this system. 
System identification is made by typing the number corre- 
sponding to the adjacent function name given in the menu 
above. In the example shown below, the Dirac delta func- 
tion is defined as the system impulse response, while the 
original signal is defined as the forcing function. 

Input transient 1 represents number ? 

>1 

Input transient 2 represents number ? 

>2 

In digital band-pass filtering, the forcing function is syn- 
onymous with the signal under investigation. The impulse 
response function represents the filter, while the transient 



response represents the filtered version of the original 
signal. By identifying two of the corresponding functions, 
as part of the linear function, the computer then calculates 
the third unknown function. In such a case, the sinusoid 
will be filtered by the Dirac delta function and output as 
the transient response function. The message below re- 
quires an indication of the domain in which to calculate 
the unknown function. 

ENTER: 1.0 (= Time-domain (de)convolution) 

2.0 (= Frequency-domain (de)convolution) 
>1 

Entering a value of 1 will cause the (de) convolution op- 
eration to take place in the time domain. Entering a value 
of 2 will utilize the frequency domain. In either case, the 
transient response is output as a function of time. The 
message appearing on the screen indicates which function 
is determined. 



***************************************************************************** 

THE TWO TIME SERIES FUNCTIONS HAVE BEEN (DE)CONVOLVED. THE 

CALCULATED FUNCTION CORRESPONDS TO NUMBER 3.00 FROM ABOVE 

***************************************************************************** 

It is this new function that is now resident in memory. 



18 



FILTER MODULE 

The use of frequency domain filters extends the flexibil- 
ity of data enhancement. By implementing the FILTER 
module (option 6), the examiner is permitted to select 
one of seven frequency-domain filters. The available 
frequency-domain filters appear in the following menu: 

FILTER MODULE 

FREQUENCY-DOMAIN FILTER MENU 



ENTER: 



1.0 | 


> BANDPASS 


with cosine window . 


) 


2.0| 


[= BANDPASS 


with half sine window 


) 


3.0< 


; = BANDPASS 


with cosine tapers . . . 


) 


4.0 I 


; = BANDPASS 


with linear tapers . . . 


) 


5.0| 


{= HIGH PASS 


with cosine tapers . . . 


) 


6.0 1 


> LOW PASS 


with cosine tapers . . . 


) 


7.0| 


{= NOTCH 


with cosine tapers . . . 


) 


8.0 1 


; = RETURN to 


menu 


) 



The investigator can employ a half-sine or cosine window, 
notch, band-pass, high-cut (low-pass), or low-cut (high- 
pass) filter. Each filter permits the user to interactively 
enter, reject, and/or pass points. With this sort of flex- 
ibility, it is possible to design the sharpness of filter cutoff. 
Both linear and cosine tapers are available. In defining 
the taper, at least two points (frequency values) must be 
defined between which the taper is implemented. In all 
cases, the taper utilized will slope from the pass point 
toward the reject point. In order for BOMSPS frequency- 
domain filters to yield realizable results, the maximum 
frequency component utilized to define a pass-or-reject (or 
reject-or-pass) combination cannot exceed the Nyquist 
frequency. To avoid a potential problem, it is suggested 
that the user examine the Fourier attribute summary table 
(option 7). Recall that the Nyquist frequency is given as 
part of the time or frequency data summary. Furthermore, 
any spectral component that is eliminated by a temporal 
filter (i.e., band-pass, low- or high-cut, etc.) operation is 
nonrecoverable. 

For a general band-pass filter, the frequency points are 
defined in the following order: low-frequency reject (lfr), 
low-frequency pass (lfp), high-frequency pass (hfp), and 
high-frequency reject (hfr). All energy outside the spe- 
cified reject points is discarded. In designing a half-sine or 
cosine band-pass filter, the user specifies only two points 
(lfr and hfr). All frequencies greater than or equal to the 
low-frequency reject point, but less than or equal to 



high-frequency reject will be passed. In defining the high- 
pass filter, the low-frequency reject or pass points are 
defined. Conversely, the low-pass filter is defined 
employing the high-frequency pass or reject points. The 
high-pass filter effectively rejects frequencies less than the 
lfr point, while the low-pass filter rejects all frequencies 
exceeding the specified hfr point. As the name "notch 
filter" suggests, all energy in excess of the lfr and less than 
the hfr is lost. Tapers are implemented automatically once 
the frequency reject or pass points are specified. 
Regardless of the type of taper (linear or cosine), energy 
defined between the low-frequency pairs and/or the high- 
frequency pairs is tapered accordingly. 

Various frequency-domain filters were applied to a 
Dirac delta function. The frequency- and time-domain 
results upon imposing these filters are displayed in fig- 
ures 12 and 13, respectively. The former figure demon- 
strates the frequency-domain equivalence of implementing 
various filters and related reject or pass points. Each 
panel is normalized and displayed from -1 to 1 along the 
abscissa, while the ordinate depicts increasing frequency 
from Hz to the folding frequency (500 Hz). While the 
actual implementation of the temporal filter occurs in the 
frequency domain, the effect is not demonstrated until an 
inverse transformation has occurred. In figure 13, the 
corresponding time-domain expressions, as filtered, are 
displayed. 

FOURIER MODULE 

When the FOURIER module (option 7) is invoked, 
BOMSPS calculates the real and corresponding imaginary 
components associated with the transformation process. 
Additionally, the amplitude, power, and phase spectra are 
determined. Phase spectral values are characteristic of 
an unwrapped sequence measured in radians. Those se- 
quence values associated with the dynamic amplitude and 
power spectra (described later) are also calculated. The 
arbitrary values are converted to a dynamic range based on 
the relation 



dB = 20 1og[s(c)|/|S( W ) max |j, (28) 



where S(w) = absolute value of transformed signal 
(spectrum). 

Consequently, dynamic range values are negative (loss of 
energy, decibels) and referenced to zero magnitude. 



19 



-10 1-10 1-10 1-10 1 



n 



(_> 



100 



200 - 



300 - 



400 - 



500 




O 

LU 
Ll_ 



100 - 



200 " 



300 ~ 



400 - 



500 




-1 1 



1 



1-10 1 



NORMALIZED AMPLITUDE 

Figure 12.-Frequency-domain filters applied to Dirac delta function (DIRAC.DAT, appendix A). Each filter is shaded to illustrate the 

portion of the spectrum that is passed. The low- and high-frequency pass and reject, in hertz, are as follows: 

Panel Low frequency High frequency 

Reject point Pass point Pass point Reject point 

A, All-pass (unfiltered spectrum) 

B, Band-pass filter with cosine window ... 500 500 

C, Band-pass filter with half-sine 

window 500 500 

D, Band-pass filter with cosine tapers .... 25 55 100 

E, Band-pass filter with linear tapers 25 55 100 

F, High-pass filter with cosine taper 25 

G, Low-frequency pass filter with 

cosine taper 55 100 

H, Notch filter with cosine tapers 25 1 00 55 



20 



0.0 



CO 



-1.00 


1.00 


A 







0.258 0.258 -0.918 0.918 -0.123 0.123, 

















i 












i 












* 












■ 












i 








1 






i 












i 






i 






i 




B 








i 







-0.133 0.133 -0.963 0.963 -0.160 0.160 -0.842 0.842 




I 
I 

l 
I 

F 





AMPLITUDE, V 

Figure 1 3.-Transient response to implementing frequency-domain filtering with regard to a Dirac delta function (DIRAC.DAT, 
appendix A). Positive polarity is shaded. The low- and high-frequency pass and reject, in hertz, are as follows: 

Panel Low frequency High frequency 

Reject point Pass point Pass point Reject point 

A, All-pass (unfiltered spectrum) 

B, Band-pass filter with cosine tapers .... 500 500 

C, Band-pass filter with half-sine 

window 500 500 

D, Band-pass filter with cosine tapers .... 25 55 1 00 

E, Band-pass filter with linear tapers 25 55 1 00 

F, High-pass filter with cosine taper 25 

G, Low-frequency pass filter with 

cosine taper 55 1 00 

H, Notch filter with cosine tapers 25 100 55 



21 



A summary table of both time and frequency attrib- is automatically displayed on the screen for the user's 
utes, associated with any waveform under investigation, information. 

TIME VARIABLES: 



SAMPLING INCREMENT 
NUMBER OF SAMPLE POINTS 
MAXIMUM TIME AMPLITUDE 



0.00100000 s 
128.00000000 
9.50000000 at 0.00400000 s 



SPECTRAL VARIABLES: 



FREQUENCY INCREMENT 
FOLDING FREQUENCY 
MAXIMUM SPECTRAL AMPLITUDE 



7.81249952 Hz 
499.99999694 Hz 
119.22966759 at 62.49999619 Hz 



*********************************************************************** 



Time variables include the sample increment, the num- 
ber of discrete sample points, and maximum relative time 
amplitude. In regard to the example given above, the 
sinusoidal function (SINUSOID.DAT) is characteristic of 
a 1-ms sample interval, 128 discrete data points, and a 
corresponding maximum amplitude of 9.5 V at 4 ms. Sim- 
ilarly, spectral variables include the frequency increment 
of 7.8 Hz, a folding frequency of 500 Hz, and a relative 



maximum spectral amplitude of 119 units at 62 Hz. The 
Fourier attributes-real transform, imaginary transform, 
amplitude, power, and phase spectra-can all be seen, in 
adjoining panels on the monitor, by affirmation of the 
following screen message. To ensure that the program is 
operating correctly, these same Fourier attributes were 
calculated and plotted for both a boxcar and a Dirac delta 
function. 



DO YOU WISH TO DISPLAY THE FOURIER ATTRIBUTES TO THE SCREEN ? 

ENTER: 1.0 (= YES) 
2.0 (= NO) 

Initially, a summary outlining the number of traces to be plotted and the folding frequency appear on the screen. 
The user is then requested to interactively set up plotting parameters. (See "Plot Module" section.) 



GAIN MODULE 

When option 8, the GAIN/VELOCITY module, is in- 
voked, the computer dispatches the following menu: 

***************** 

GAIN MODULE 

***************** 



ENTER: 1.0 (= Convert voltage gain (dB) ) 

2.0 ( = Convert voltage to velocity values . . ) 

3.0 (= Return to menu ) 

>1 

ENTER: Gain (in decibels, dB) = ? 



The GAIN/VELOCITY module has twofold impor- 
tance as a signal processing tool. First, it is possible to 
increase or decrease the gain (decibels). A typical use for 
this operation is to recover the true gain of a signal 
recorded in the field. Second, the GAIN/VELOCITY 
module can be utilized for converting from voltage to 
velocity values (or velocity to voltage values). In the ex- 
ample given above, a gain value of 1 dB is introduced. 
After each operation is finished, the menu reappears on 
the screen. The second application (amplitude conversion) 
can be implemented by selecting 2. Immediately following 
this selection, the user is requested to define the type of 
transducer response function (volts per inch per second) to 
be used. 



>50 



>2 

DO YOU WISH TO CONVERT VOLTAGES TO VELOCITY VALUES ? 



ENTER: 1.0 (= YES) 
2.0 (= NO) 
>1 



22 



ENTER: Transducer response function (vlts/in/sec) = ? 

1.0 (= single value) 

2.0 (= external function) 
>1 

ENTER: Scalar response = ? 

>1.5 

If the single-value option is selected, then one value is 
entered. Such a case is useful if the signal spectrum oc- 
curs within the transducer's flat response spectrum. How- 
ever, if the broadband energy also occurs outside the 
transducer's flat response, the investigator will be more 
likely to choose the second option (input a response func- 
tion). If the second option is selected, the user must enter 
the name of the appropriate external file that holds the 



complete transducer function. This transducer response 
function must be compatible with the signal under 
investigation. 

HILBERT MODULE 

Provisions in BOMSPS permit the calculation of trace 
attributes such as the envelope (instantaneous amplitude), 
quadrature (Hilbert transform), instantaneous phase, in- 
stantaneous frequency, and derivative. Although these 
functions represent time-domain attributes, each calcula- 
tion is simplified by conducting operations in the frequency 
domain. Upon selection of the HILBERT module (option 
9), each of these trace attributes is calculated. Immedi- 
ately following the selection of the Hilbert attribute option, 
a message will appear on the screen inquiring as to wheth- 
er these functions should be plotted on the screen. 



DO YOU WISH TO DISPLAY THE HILBERT ATTRIBUTES TO THE SCREEN ? 

ENTER: 1.0 (= YES) 
2.0 (= NO) 
>1 



As is the case for the Fourier attributes, the trace attrib- 
utes are displayed in adjoining panels after the plot setup 
procedure is completed. Figure 14 depicts the screen 
view of those Hilbert attributes related to the sinusoid 
(SINUSOID.DAT, appendix A). In this example, the 
100 data points are each sampled at 1-ms intervals. The 
first panel depicts the resident waveform (input signal). 
The remaining panels (from left to right) are the quad- 
rature, envelope (instantaneous amplitude), and instan- 
taneous phase. Time-amplitude values are plotted with 
time along the ordinate and amplitude along the abscissa. 
Again, the abscissa values are indicative of a normalized 
set. 

Basic operations such as gain control, muting, normal- 
ization, polarity change, randomization, scaling, shifting 
(translation), stacking, smoothing, and windowing account 
for the majority of time-domain digital enhancement (fil- 
tering) operations employed by BOMSPS. Routine use 
of these operations, prior to implementation of other fil- 
ters and/or display features, warrants their inclusion. The 
utilization of any of these menu operations results in re- 
placing the original contents (in memory) with a perturbed 
(filtered) version. If, at any time, the user wishes to recall 
the original time series function into memory, he or she 
simply invokes the reset operation (option 15). This reset 
process results in reloading the temporary memory with 
the original waveform held in a permanent memory reg- 
ister. Upon completion of this process, the following 
message appears on the screen indicating such: 

*********************************************** 

RESET MEMORY: Original Time Series Function 

*********************************************** 



Higher order time-domain processing involving convolution 
and/or cross correlation is useful but requires an operator 
(second time function) in order to be executed. Both of 
these operations (convolution and cross correlation) are 
discussed in detail under the heading "Two Time Series 
Functions." 

INTERPOLATE MODULE 

For data to be processed utilizing BOMSPS, discrete 
values must be sampled at equal time intervals. If a 
graphics tablet is employed, the signal generally ends up 
digitized in a random fashion with unevenly spaced time 
values. Any attempt to employ options in the menu other 
than option 10 will result in an error message and termi- 
nation of the program. As an example, consider reeval- 
uating a data file characteristic of even sample intervals, 
but having too few data points (86). The dispatched mes- 
sage is 



>10 



******************* 

INTERP MODULE: 

******************* 

ENTER: Number of Points (numpts): Total number of 
new values to be interpolated and stored. 

>1000 



23 



0.£0 



Input 3ignal 



Quadrature 



Instant. Amplitude Instant. Pha3e 



v> 



UJ 




-1 



-1 



-1 



-1 



AMPLITUDE 

Figure 1 4.-Screen view of trace attributes: original transient (input signal), quadrature, instantaneous amplitude 
(envelope), and instantaneous phase functions. 



In the example, the examiner chooses to increase the 
sample density from 86 data points (sampled at 1-ms in- 
tervals) to 1,000 data points (sampled at 86-/xs intervals). 
The previous number of 86 discrete data pairs (time, am- 
plitude) is now replaced by 1,000. Immediately following 
this entry is a message regarding the uniformity of sam- 
pling (below). 

TYPE: 1.0 (= YES) 
2.0 (= NO) 
>1 



The default response (no) implies that the data set is of 
random origin. The user would then enter the required 
sample interval (seconds) to be employed during the in- 
terpolation. In the example given, however, it is of interest 
to reevaluate an evenly spaced data set. Either a spline or 
quasi-spline routine (5) can be employed for this analysis. 
Depending on the complexity of a given waveform, one 
interpolation routine may yield better results than another. 
Hence, there exists the availability of two different in- 
terpolation routines. The example given uses the cubic 
spline. 



PLEASE SELECT ONE OF THE FOLLOWING OPERATIONS 



TYPE: 1.0 (= CUBIC SPLINE) 

2.0 (= QUASI-HERMITE SPLINE) 
>1 



24 



A cubic spline is implemented by entering a value of 1. 
Regardless of the user's selection, both interpolation rou- 
tines result in returning an evenly spaced set of discrete 

ERROR MESSAGE 



data points. Following the appropriate selection, BOMSPS 
returns a table of information. 



x(i) y(i) c(i,l) c(i,2) c(i,3) 



1 


0.0000 


0.0000 


-5756.7339 


0.6920 


0.0000 


2 


0.0010 


0.0000 


4303.3667 


0.6930 


0.0000 


3 


0.0020 


5.7000 


5643.2681 


0.6940 


0.0000 


4 


0.0030 


9.5000 


1623.5597 


0.6950 


0.0000 


5 


0.0040 


9.1200 


-1877.5096 


0.6960 


0.0000 


6 


0.0050 


6.2300 


-3923.5193 


0.6970 


0.0000 


7 


0.0060 


1.5900 


-5018.4102 


0.6980 


0.0000 


8 


0.0070 


-3.2400 


-4412.8354 


0.6990 


0.0000 


9 


0.0080 


-6.8600 


-2680.2505 


0.7000 


0.0000 


10 


0.0090 


-8.4300 


-436.1656 


0.7010 


0.0000 



DO YOU WISH TO SEE THE REMAINING DATA PTS ? 

TYPE: 1.0 (= YES) 
2.0 (= NO) 
>2 



The first line denotes an error message number. If a 
message other that zero appears, there is some problem in 
the data file. A value of 129 denotes that the row dimen- 
sion is less than the total number of discrete points. If the 
number 130 appears, then the number of elements in the 
external data file is less than the new number. An error 
message of 131 indicates that input abscissa values are not 
ordered in a successive manner. Directly below the error 
message are written six columns of information. The first 
column (i) depicts the relative position of a data pair. The 
second and third columns are the time x(i) and relative 
amplitude y(i), respectively. The last three columns (cl, 
c2, c3) represent matrix coefficients used in the interpo- 
lation polynomials. The first 10 lines of information are 
automatically written to the screen, for quality control. 
However, the user may opt to see the remaining data by 
entering a 1. The trailing information indicates where the 
reevaluation data can be found. This data file is formatted 
such that it can be read directly by BOMSPS when the 
following file name is typed in: 

> INTERP1.DAT 

Other available information includes the total time dura- 
tion of the signal, the original number of discrete data 
pairs, the number of reevaluated data pairs, and the num- 
ber of reevaluated data pairs padded to the next power of 
2. The power-of-2 criterion is necessary for transforming 
data. A more complete explanation is given in the section 
"Numerical Considerations." A more complete explanation 



of these interpolation routines can be found in the doc- 
umentation provided by International Mathematics and 
Statistical Library (IMSL) (5). Since both of these 
interpolation routines are called from the IMSL, the user 
must have this available for the operation to be imple- 
mented successfully. If the IMSL library is not available, 
the call line at which BOMSPS employs INTERP must be 
commented out and recompiled. 

Figure 15 represents the results of implementing the 
interpolation operation (INTERP module, option 10). 
Figure 15A presents a discrete version of the original time 
series function, composed of 100 data points sampled at 
1-ms intervals. Despite the increased sample density, the 
signal maintains the same shape in figure 15B. 

MUTE MODULE 

Another commonly utilized operation is the muting of 
a waveform. Selecting the MUTE module (option 11) 
results in dispatching the following message: 

MUTE MODULE 

ENTER: Beginning time for mute window = ? 
>0 

ENTER: End time for mute window = ? 
> 0.030 



25 



•9.4 9.4 -9.4 9.4 -1.0 



Q 



Cl 
H 
< 

LU 
> 

< 

_1 
LU 





□ 




1 1 1 

KEY A 




D 

□a 




a Sample point 


.4 


- □ 

a a 






.2 


D 




r\ * - 




i_ _u 




D g D ° D 3 










-.2 


° c 
□ 


a 

D D 


V 




□ 


□ ° 




-.6 


a 
□ 

D 


erf 





1.0 



-2 



1.0 



1 


1 1 

B 


oB A 




a° W* 




*° fl 


- 


D D R H 




n D n Pi — 








a a 5 « # 




a a B 1 / 


«. V% 


n S S 1 f 




a a a fl fl 




3- — fl — Q — §— — — 

s ni r 


"^ ~ — m ~\ ^^^^^^^ 


- S a 1 # 


\ § W 


Ss / 




° B I / 




^1 l/ 




1 I 





0.02 



0.04 



0.06 



0.08 



0.1 



TIME, s 

Figure 15.-lnterpolation of a digital time sequence 
(SINUSOID.DAT, appendix A). A, Original signal, 100 discrete 
data points sampled at 1-ms intervals; B, interpolated signal, 
1,000 discrete data points sampled at 100-ms intervals. 



The user is required to enter first the beginning and 
ending time points of the mute window (seconds). The 
window chosen in the example is 30 ms. This effectively 
zeros all time points within this mute window (including 
the endpoints). The result of using a 30-ms mute window 
on the original sinusoid shown in figure 16A is illustrated 
in figure 16B. This operation can be used effectively when 
the actual analysis might involve only a portion of a com- 
plete wave train. E.g., upon collecting a complete seismic 
wave train, the investigator may wish to analyze only the 
surface wave (removing the compressional and shear 
components). 



.05 



.05 - 



.05 - 




•9.4 



9 4 



1.0 



AMPLITUDE, V 



Figure 16.-Some time-domain operations as applied to a time 
series function (SINUSOID.DAT, appendix A). A, Input signal; B, 
mute window: from to 30 ms; C, normalized; D, 180° phase 
reversal; E, randomized; F, time shift: 30 ms, G, five-point 
smooth applied to 15E; H, cosine window: from to 100 ms. 

NORMALIZATION MODULE 

In some instances, it may be important to normalize 
the results. This is easily handled by invoking the nor- 
malization operation (NORM module, option 12). The 
NORM module menu (below) permits the user to normal- 
ize any or all of the functions calculated by dividing the 
entire array of values by the maximum amplitude. 

****************** 

NORM MODULE 

****************** 



ENTER: 



>2 



1.0 (= Normalize all functions) 
2.0 (= Normalize some functions) 
3.0 (= Return to menu) 



26 



All functions can be normalized by entering a value of 1. 
If it is of interest to normalize more than one function but 
less than the total number, the user simply selects the 
second option. This action results in a scrolling through a 
series of messages inquiring as to whether the function 
should be normalized. The example below depicts the 
series of interactive questions pertaining to selective nor- 
malization of the functions. A normalized version of the 
time series transient (SINUSOID.DAT, appendix A) is 
included in figure 16C. 

NORMALIZE: TRANSIENT (REAL) ? 

ENTER: 1.0 (= YES) 
2.0 ( = NO) 
>1 

NORMALIZE: INPUT TRANSIENT (IMAG.) ? 
>2 

NORMALIZE: CUMULATIVE ENERGY (TIME) ? 

>2 



PLOT: Every Nth data point; N = ? 
>1 

MINIMUM: Y-axis cutoff point = ? 
ENTER: (= NONE) 

N (= ARBITRARY POINT) 
>0 

MAXIMUM: Y-axis cutoff point - ? 
ENTER: (= NONE) 

N (= ARBITRARY POINT) 
>500 

INCREMENT: Y-axis time/freq= ? 
ENTER: (= NONE) 

N (= ARBITRARY POINT) 
>100 

GRID LINES: Y-axis time/freq = ? 
ENTER: (= NONE) 

1 (= YES) 
>1 



NORMALIZE: ENVELOPE (INST. AMPLITUDE) ? 
>1 

NORMALIZE: QUADRATURE ? 

>2 

NORMALIZE INSTANTANEOUS PHASE ? 

>2 

NORMALIZE: TRANSFORM SPECTRUM (REAL) ? 
>1 

NORMALIZE: TRANSFORM SPECTRUM (IMAG.) ? 
>1 

NORMALIZE: AMPLITUDE SPECTRUM ? 
>1 

NORMALIZE: PHASE SPECTRUM ? 
>1 

PLOT MODULE 

The prompts dispatched to the screen are depicted 
below. This information is useful for establishing the 
minimum or maximum display coordinates and other rel- 
evant display parameters. 

>1 

PLOT MODULE: 

NUMBER OF TRACES TO BE PLOTTED = 4 

FOLDING FREQUENCY = 500.0 

********************************************* 



VARIABLE AREA: Shading = ? 
ENTER: (= NONE) 

1 ( = YES) 
>1 

The plotting speed is governed primarily by the baud rate, 
the bit map process (as opposed to a raster mapping pro- 
cess), and the relative number of data points in the array 
under investigation. For these reasons, it is recommended 
that the baud rate be set to the fastest possible (i.e., 9,600 
for most cases). When many data points exist, the user 
may increase the relative plotting speed by choosing to plot 
every Nth data point. In the example given, every data 
point will be displayed on the screen. The next two mes- 
sages refer to setting the window minimum and maximum 
time or frequency values. Here the minimum and maxi- 
mum frequency values correspond to and 500 Hz, re- 
spectively. The upper limit (maximum frequency) is es- 
tablished by the folding frequency given in the summary. 
The frequency increment in this case is set equal to 
100 Hz. To further enhance the plot, grid lines (drawn at 
every increment) and/or shading may be implemented. 
The shading (variable area) occurs above the point where 
the transient quantity crosses the zero axis. If time or 
frequency lines are not chosen, the program defaults to 
tick marks written along the ordinate. In the example, 
both the time or frequency lines and shading are chosen. 
The setup of display parameters permits flexibility in an- 
alyzing a small window, if need be. Figure 17 depicts the 
Fourier attributes characteristic of an arbitrary time series 
function (SINUSOID.DAT, appendix A). The display 
character reflects those setup parameters chosen above. 
The ordinate values are frequencies (from to 500 Hz). 
Each corresponding function is normalized, with the pos- 
itive values shaded in. Each spectrum panel is labeled and 
displays its normalized equivalent. Detailed information 



27 



Rea I ( transform) I mag < transform ) 



Rmp I i tude 



Phase 



N 

u 

-Z- 
LU 

Z) 

o 




100 



200 _. 



300 



400 



500 



-1 







-1 
1 



-1 







1 



AMPLITUDE 

Figure 17.-Screen view of Fourier attributes; transformation (real and imaginary), amplitude, and phase spectra. 



relative to a given function amplitude can be ascertained 
by invoking the display operation (option 4 in the main 
menu). It should be noted that data can also be written to 
scratch files in ASCI format (option 22) and plotted using 
a commercially available graphing package. 

POLARITY MODULE 

Invoking the POLARITY module (option 13) results in 
the following menu: 

POLARITY MODULE 

********************** 

ENTER: 1.0 (= 180 degree phase modification) 
2.0 ( = 90 degree phase modification) 
3.0 (= Return to menu ) 

>1 

In this case, the user can modify the signal's phase by 
either 90° or 180°. Figure 16D depicts a 180° phase re- 
versal of a signal (SINUSOID.DAT, appendix A). A 



practical use of this option is in determining the first break 
associated with the shear wave arrival. In such a case, two 
shear waves are superimposed (one normal and one 180° 
out of phase). From such an analysis, it is possible to 
more accurately assess the compressional wave's travel 
time and, indirectly, the velocity. In the interpretation of 
a two-dimensional seismogram, reflection continuity can 
sometimes be favorably enhanced upon converting the 
polarity. 

RANDOMIZATION MODULE 

For the sake of instructional purposes, BOMSPS in- 
corporates a random number generator (RANDOM mod- 
ule, option 14). If the randomization operation is utilized, 
a message will appear on the screen when it is completed. 

****************************** 
RANDOM MODULE: finished 



Typically, the random number generator is used to ex- 
amine the usefulness of various band-pass filters. At a 



28 



later point in the text, the randomized signal is used to 
verify the usefulness of the smoothing operation. Fig- 
ure 16E depicts the randomized version of the sinusoid. 

SCALE MODULE 

The SCALE module (option 16) permits an investigator 
to interactively choose to scale either the abscissa or 

Hi***************** 
SCALE MODULE 



ordinate values of a given input signal. Scaling occurs by 
either adding a scalar constant to each amplitude or time 
value, or multiplying each value by a constant. In the 
example given below, the amplitude values (Y-axis) are 
multiplied by a factor of 1. Since it is decided not to scale 
the time (X-axis) values, the signal processing menu ap- 
pears on the screen. 



DO YOU WISH TO SCALE THE AMPLITUDE (Y) VALUES 

ENTER: 1.0 (= YES) 
2.0 (= NO) 
>1 

ENTER: Scaling constant = ? 
>1 

ENTER: 1.0 (= additive) 

2.0 (= multiplicative) 

>2 

DO YOU WISH TO SCALE THE TIME (X) VALUES ? 

ENTER: 1.0 (= YES) 
2.0 (= NO) 
>2 



SHIFT MODULE 

Shifting (or translation) of a sampled data series can be 
facilitated by utilizing either the SHIFT (or TRANSLATE) 
module found in the BOMSUBS library (option 17). The 
former module operates in the frequency domain, by mod- 
ifying the phase. The latter module operates in time, 
perturbing the time values by a scalar constant, much the 
same as the last option. In the main program, a call line 
is incorporated for each of these subroutines. Since only 
one of these routines can be used, the other routine is 
commented upon. The user may implement the alternate 
module by swapping comments at the call sequence in the 
main program. The program must then be recompiled. 
Upon selecting the SHIFT operation (option 17), the user 
enters the amount of time (in seconds) that the trace is to 
be translated. In the example given below, a shift of 20 ms 
is entered. 

SHIFT MODULE 



ENTER: 



>1 



1.0 (= TIME-SHIFT TRANSIENT) 
2.0 (= RETURN TO MENU . . . ) 



ENTER: Time shift duration ( + /- sec) = ? 
> 0.020 

The effect of imposing a 20-ms (0.020-s) time shift is 
depicted in figure 16F. 

In either time b(t-to) or frequency B(o>-wo), a shift 
results in multiplying the inverse Fourier transform of the 
unshifted version by an exponential term. The time-shifted 
Fourier transform pair is given by 



B(w) = exp (jwt +/- wT ). 



(29) 



Since time shifting of a function does not alter the 
magnitude of the Fourier transform, it immediately follows 
that the relative amplitude spectrum would also remain 
invariant to a time shift. The effect of time shifting, 
however, does modify the phase spectrum 4>(u). Hence, 
the transform (or amplitude) spectrum is nonunique and 
is not sufficient information to define the input time 
function. In other words, two signals differing only in 
phase (nonunique) spectra will be unrelated in time. This 
fact is true despite their similarity in amplitude spectra. 
The original phase wt is effectively modified by adding a 
linear ramp wTo. By adding a linear ramp (function of 
frequency) to the phase spectrum and taking the inverse 
Fourier transform, the original signal will be translated 
along the abscissa (time axis) by To. 



29 



Two time functions that have the same phase spectra, 
but whose zero-frequency terms differ, will appear only to 
be shifted (positive or negative) along the ordinate. As 
translation increases along the ordinate (parallel to the 
abscissa), the time function becomes physically suppressed. 
Eventually, the signal approaches that of a boxcar function. 
This phenomenon is known as the time-domain windowing 
effect and must be considered when scaling a transient. 
The effect is to suppress the higher frequency energy, 
ultimately approaching a Sine function. To reduce the 
unwanted suppression of the original transient, the time 
series should be padded with additional data points (zeros) 
as illustrated in figure 18. As a rule of thumb, the total 
number of discrete data points including zeros should 
represent twice the next power of 2. Additionally, one 
may wish to taper the signal, in order to promote side lobe 
suppression. Tapering can be done by implementing any 
one of the available frequency-domain filters discussed 
under the heading "Filter Module." 



SMOOTH MODULE 

A particularly useful operation for data enhancement, 
prior to display, is the smoothing operation (option 19). 
The SMOOTH module provides either a three- or a five- 
point running window, over which values are averaged. 
The mathematical expressions used for these operations 
are- 
Three-point smooth: 



i + 2 



s(t) i + 1 = Ss(t)/3i = l,2,...,(N-l). 
j = l 



(30) 



Five-point smooth: 



s(t) i+2 =!^s(t)/5 i = l,2,...,(N-3). 



(31) 



1 ■ 



Q_ 


2 


z: 


< 


1 


LU 




> 





1— 




< 




_J 




LU 


3 


QL 




2 




1 








;, aa 



i - 



-WV^ 



^W 



TIME 



Figure 18. -Time-domain windowing cause, effect, and solution. A, Originally sampled time sequence; 
B, translation along ordinate (y-axis). Note original shape approaches a boxcar function. C, Padding the 
signal to next higher power of 2, e.g., 2 s to 2* data points, reduced the amplitude of the boxcar. 



30 



The appropriate smoothing operation can be initiated by 
entering a value of 1 (three-point) or 2 (five-point). In 
the example given, a three-point smoothing operation is 
initiated. 



SMOOTH MODULE 



ENTER: 1.0 (= 3 point smoothing) 

2.0 ( = 5 point smoothing) 

3.0 (= Return to menu. .) 
>1 

The result of applying this smoothing operation to the 
randomized transient (fig. 16E) is depicted in figure 16G. 

STACK MODULE 

In many instances, the amplitude of an individual time 
transient may be readily discernible above background 
noise. The stacking process (STACK module, option 18) 
results in adding successive time functions (multiple rec- 
ords). More exactly, a time point of the original signal is 
added to the corresponding time point of the successive 
waveform. This process improves the signal-to-noise ra- 
tio. The investigator should be aware, however, that the 
relative frequency content will decrease. In the example 
outlined below, two signals, DIRAC.DAT and 
BOXCAR.DAT, are added to the original transient. 

****************** 

STACK MODULE 

****************** 



ENTER: Number of additional waveforms to 
be stacked with the original = ? 

>2 

ENTER: External filename = ? 

> dirac.dat 

ENTER: External filename = ? 

> boxcar.dat 

WINDOW MODULE 

********************* 



WINDOW MODULE 

Windowing (shaping) of a given waveform is used con- 
ventionally as a filter scheme, with the intention of sup- 
pressing side lobe energy. Windowing of the time function 
may be conducted in BOMSPS by implementing either a 
half-sine, Hanning, or cosine weighting function. These 
weighting functions are incorporated into the WINDOW 
module and can be invoked by selecting option 20. De- 
pending on which window is chosen, the resulting function 
is based on one of the following relationships: 



The half-sine weighting function W hs (t) is 

W hs (t) = sin(t/T); 
the cosine weighting function W c (t) is 



W c (t) = 0.5 



l-cos(2t/T) 



(32) 



(33) 



the Hanning weighting function W h (t) is 



W h (t) = 0.5 



l + cos(t/T) 



(34) 



where t = time 



and 



T = period of window. 



The half-sine, cosine, and Hanning window's theoretical 
rolloffs are 12, 18, and 24 dB per octave, respectively. This 
is a useful guide in determining the appropriate window to 
use. While shaping tends to increase the signal-to-noise 
ratio, it does so at the expense of higher relative frequen- 
cies. More specifically, shaping tends to remove relatively 
high frequency energy. 

When the WINDOW module is invoked, a menu ap- 
pears listing the various weighting functions available. In 
the example given below, a 100-ms Hanning window is 
used. The effect of utilizing this window on the original 
transient is depicted in figure 16H. 



ENTER: ********** WEIGHTING FUNCTION *********** 

1. Hanning function (18 dB/octave) 

2. Half-sine function (24 dB/octave) 

3. Cosine function (30 dB/octave) 

4. Return to menu ) 

>1 

ENTER: Beginning time point (sec) = ? 
>0 

ENTER: End time point (sec) = ? 
>0.10 



WRITE MODULE 

It is often of interest to examine individual data points 
prior (or in addition) to displaying them graphically. Cal- 
culated time- and/or frequency-domain functions can be 



31 



displayed on the screen by invoking the WRITE operation 
(option 21) from the signal processing menu. Below is the 
first of two questions, prompting the user as to whether to 
write the time-domain data to the screen. 



WRITE TIME-DOMAIN DATA TO THE SCREEN ? 

ENTER: 1.0 (= YES) 
2.0 (= NO) 
>1 



An affirmative response results in the following summary table: 
TIME DOMAIN ATTRIBUTES: 



TIME 


SIGNAL 


CUM.ENERGY 


QUADRATURE 


ENVELOPE 


INST.PHASE 


0.00000 


0.00000 


0.00000 


0.00000 


0.00000 


0.00000 


0.00100 


0.00000 


0.00000 


0.00000 


0.00000 


0.00000 


0.00200 


0.05048 


0.00000 


0.00000 


0.00000 


0.00000 


0.00300 


0.14923 


0.00000 


0.00000 


0.00000 


0.00000 


0.00400 


0.22318 


0.00000 


0.00000 


0.00000 


0.00000 


0.00500 


0.21875 


0.00000 


0.00000 


0.00000 


0.00000 


0.00600 


0.07566 


0.00000 


0.00000 


0.00000 


0.00000 


0.00700 


-0.20038 


0.00000 


0.00000 


0.00000 


0.00000 


0.00800 


-0.53395 


0.00000 


0.00000 


0.00000 


0.00000 


0.00900 


-0.80499 


0.00000 


0.00000 


0.00000 


0.00000 



Note that functions to the right of the signal (cumulative energy, quadrature function, envelope, and instantaneous phase) 
are assigned values of zero. This will occur if the Hilbert attribute operation was not previously selected. The next 
question requires the user to determine whether the frequency-domain data should be written to the screen. 

WRITE FREQUENCY DOMAIN DATA TO SCREEN ? 
ENTER: 1.0 ( = YES) 
2.0 (-NO) 

In the example, the first of 10 frequency-domain attributes was investigated. The table of these Fourier attributes is 
given below. 

FOURIER-DOMAIN ATTRIBUTES: 

FREQUENCY TRANSFORM AMPLITUDE POWER PHASE 



0.00 


61.32 


61.32 


3760.14 


0.00 


7.81 


-19.69 


47.55 


2271.90 


2.00 


15.62 


7.06 


71.18 


5066.84 


-1.47 


23.44 


90.98 


99.49 


9899.19 


0.42 


31.25 


-49.23 


73.08 


5340.83 


2.31 


39.06 


-2.37 


39.95 


1596.16 


-1.63 


46.87 


6.75 


51.25 


2627.33 


-1.44 


54.69 


119.22 


119.23 


14215.78 


-0.01 


62.50 


55.34 


73.90 


5461.31 


0.72 


70.31 


34.86 


112.03 


12550.75 


1.25 



In the above summary table, all of the functions listed display spectral information. If, however, the Fourier attribute 
option is not selected prior to writing this information to the screen, all values would be zero. 



32 



Provided that the data are satisfactory, the user may 
transfer any or all of the functions calculated to external 
fdes. Prior to data transfer and storage to external files, 
two conditions must be met. First, the user must ensure 
that the correct operators are selected. These WRITE 
operators (1 = write, 2 = do not write) are entered in the 
fde WRT.IN, which in turn is read by BOMSPS (see ap- 
pendix A). Proper selection of WRITE operators permits 



their duplication to external files, only if the second" cri- 
terion is met. To initiate data transfer, the user must 
select option 22 (WRITE). The resulting menu from the 
selection of option 22 permits the user to interactively 
determine the total number of data pairs to be written to 
the scratch files. In the example below, the total number 
of discrete points is written to each external data fde. 



HOW MANY PAIRS OF DISCRETE DATA POINTS ARE TO 
BE WRITTEN TO A SCRATCH FILE(S) ? 

ENTER: 1.0 (= Arbitrary number of data points) 

2.0 (= Total number of data points (NP2)) 
>1 



ENTER: 

>128 



Number of points = ? 



Any one of these functions can be graphically displayed on 
the screen (prior to writing out data to a fde) by invoking 
the DISPLAY module (option 4). The same setup menu, 
employed for both Fourier and Hilbert attributes, appears 
on the screen. The primary difference, as compared with 
the former displays of Fourier and Hilbert attributes, is 



that only one function can be displayed. However, addi- 
tional detad is included with respect to amplitude values 
(unnormalized) along the abscissa. In the former displays, 
the amplitude is normalized (unity). For this reason, the 
amplitude is not displayed along the abscissa. 



SUMMARY 



The Bureau has advanced the state-of-the-art signal 
processing for microcomputers with the development of 
BOMSPS. This information circular detads those signal 
processing concepts and expressions incorporated into 
the BOMSPS code and provides a user's document for its 
implementation. The Fourier transformation and its 
properties are reviewed to dlustrate applicability in 



developing those mathematical expressions used as digital 
operators. The advantages and limitations of basic and 
higher order time- and frequency-domain operations are 
discussed. The examples and easy-to-follow menus out- 
lined in the tutorial facilitate the implementation of 
BOMSPS on a routine basis. 



REFERENCES 



1. Brigham, E. The Fast Fourier Transform. Prentice-Hall, 1974, 
248 pp. 

2. Bracewell, R The Fourier Transform and Its Applications. 
McGraw-Hill, 1965, 383 pp. 

3. Champeneny, D. Fourier Transforms and Their Physical 
Applications. Academic, 1973, 265 pp. 



4. Robinson, E. Statistical Communication and Detection With 
Special Reference to Digital Data Processing of Radar and Seismic 
Signals. Chades Griffin Co., London, 1967, 362 pp. 

5. International Mathematics and Statistics Library (Houston, TX). 
Reference Manual. V. 3, ch. 1, 1982. 



33 

APPENDIX A.-SAMPLE DATA FILES 

Sample data files are shown for the following functions: and relative amplitude values— is indicated for the first data 

boxcar (BOXCAR.DAT), Dirac delta (DIRAC.DAT), and file. Also, a sample data file for the write operator 

sinusoid (SINUSOID.DAT). The format-including the (WRT.IN) is included, 
total number of discrete data pairs, sample interval, time, 

Filename: BOXCAR.DAT 



(DELT = sample 


(NP1 = total number of 


interval) 


discrete data pairs) 


.001 


86.00 


(DLT(NPl) = 


(SI (NP1) = relative 


time (seconds)) 


amplitude) 


.000 


1.00 


.001 


1.00 


.002 


1.00 


.003 


1.00 


.004 


1.00 


.005 


1.00 


.006 


1.00 


.007 


1.00 


.008 


1.00 


.009 


1.00 


.010 


0.00 


.011 


0.00 


.012 


0.00 


.013 


0.00 


.014 


0.00 


.015 


0.00 


.016 


0.00 


.017 


0.00 


.018 


0.00 


.019 


0.00 


.020 


0.00 


.021 


0.00 


.022 


0.00 


.023 


0.00 


.024 


0.00 


.025 


0.00 


.026 


0.00 


.027 


0.00 


.028 


0.00 


.029 


0.00 


.030 


0.00 


.031 


0.00 


.032 


0.00 


.033 


0.00 


.034 


0.00 


.035 


0.00 


.036 


0.00 


.037 


0.00 


.038 


0.00 


.039 


0.00 


.040 


0.00 


.041 


0.00 


.042 


0.00 



34 



File name: BOXCAR.DAT~Continued 

(DLT (NP1) = (SI (NP1) = relative 

time (seconds)) amplitude) 

.043 0.00 

.044 0.00 

.045 0.00 

.046 0.00 

.047 0.00 

.048 0.00 

.049 0.00 

.050 0.00 

.051 0.00 

.052 0.00 

.053 0.00 

.054 0.00 

.055 0.00 

.056 0.00 

.057 0.00 

.058 0.00 

.059 0.00 

.060 0.00 

.061 0.00 

.062 0.00 

.063 0.00 

.064 0.00 

.065 0.00 

.066 0.00 

.067 0.00 

.068 0.00 

.069 0.00 

.070 0.00 

.071 0.00 

.072 0.00 

.073 0.00 

.074 0.00 

.075 0.00 

.076 0.00 

.077 0.00 

.078 0.00 

.079 0.00 

.080 0.00 

.081 0.00 

.082 0.00 

.083 0.00 

.084 0.00 

.085 0.00 



35 



Filename: DIRAC.DAT 

.001 86.00 

.000 1.00 

.001 0.00 

.002 0.00 

.003 0.00 

.004 0.00 

.005 0.00 

.006 0.00 

.007 0.00 

.008 0.00 

.009 0.00 

.010 0.00 

.011 0.00 

.012 0.00 

.013 0.00 

.014 0.00 

.015 0.00 

.016 0.00 

.017 0.00 

.018 0.00 

.019 0.00 

.020 0.00 

.021 0.00 

.022 0.00 

.023 0.00 

.024 0.00 

.025 0.00 

.026 0.00 

.027 0.00 

.028 0.00 

.029 0.00 

.030 0.00 

.031 0.00 

.032 0.00 

.033 0.00 

.034 0.00 

.035 0.00 

.036 0.00 

.037 0.00 

.038 0.00 

.039 0.00 

.040 0.00 

.041 0.00 

.042 0.00 

.043 0.00 

.044 0.00 

.045 0.00 

.046 0.00 

.047 0.00 

.048 0.00 

.049 0.00 

.050 0.00 

.051 0.00 

.052 0.00 

.053 0.00 

.054 0.00 

.055 0.00 



36 



Filename: DIRAC.DAT— Continued 

.056 0.00 

.057 0.00 

.058 0.00 

.059 0.00 

.060 0.00 

.061 0.00 

.062 0.00 

.063 0.00 

.064 0.00 

.065 0.00 

.066 0.00 

.067 0.00 

.068 0.00 

.069 0.00 

.070 0.00 

.071 0.00 

.072 0.00 

.073 0.00 

.074 0.00 

.075 0.00 

.076 0.00 

.077 0.00 

.078 0.00 

.079 0.00 

.080 0.00 

.081 0.00 

.082 0.00 

.083 0.00 

.084 0.00 

.085 0.00 

File name: SINUSOID.DAT 

.001 86.00 

.000 0.00 

.001 0.00 

.002 5.70 

.003 9.50 

.004 9.12 

.005 6.23 

.006 1.59 

.007 -3.24 

.008 -6.86 

.009 -8.43 

.010 -7.78 

.011 -5.31 

.012 -1.80 

.013 1.86 

.014 4.90 

.015 6.82 

.016 7.39 

.017 6.70 

.018 5.03 

.019 2.74 

.020 .26 

.021 -2.07 



37 



File name: SINUS OID. DAT-Continued 

•022 -3.98 

.023 -5.32 

.024 -6.04 

.025 -6.17 

.026 -5.79 

.027 -5.01 

.028 -3.98 

.029 -2.81 

.030 -1.59 

.031 -.43 

.032 .63 

.033 1.56 

.034 2.33 

.035 2.95 

.036 3.43 

.037 3.77 

.038 4.01 

.039 4.15 

.040 4.23 

.041 4.25 

.042 4.24 

.043 4.19 

.044 4.13 

.045 4.06 

.046 3.98 

•047 3.89 

•048 3.79 

•049 3.68 

-050 3.54 

•051 3.38 

.052 3.18 

.053 2.94 

.054 2.65 

.055 2.30 

.056 1.89 

.057 1.42 

.058 .88 

•059 .30 

.060 -.31 

•061 -.92 

.062 -1.51 

•063 -2.02 

.064 -2.42 

.065 -2.42 

•066 -2.70 

.067 -2.52 

.068 -2.10 

•069 -1.47 

.070 -.68 

•071 .19 

•072 1.04 

•073 1.73 

•074 2.17 

•075 2.25 

•076 1.95 

•077 1.30 

.078 .40 



38 



File name: SINUSOID. DAT-Continued 

.079 -.57 

.080 -1.41 

.081 -1.90 

.082 -1.92 

.083 -1.44 

.084 -.58 

.085 .44 

Filename: WRT.IN 

2.0002.0002.0002.0001.0002.0002.0002.000 
2.0002.0001.0002.0002.0002.0002.0002.000 



39 

APPENDIX B.-MATHEMATICAL RELATIONSHIP OF FOURIER 
TRANSFORMATION TO HILBERT TRANSFORMATION 

This relationship is used to develop the Hilbert attributes (i.e., quadrature, instantaneous amplitude, and instantaneous 
phase functions). 



By definition, 



A A 

s(t) +♦ -j Sgn(w) S(u>) = S(w), 



A 

where s(t) = Hilbert transform of a time series function s(t) (quadrature function), 

t = time, 

** = Fourier transform pair, 

j = (-1)*, 

Sgn(w) = Signum function, 

w = angular frequency, 

S(w) = Fourier transform of s(t), 

A A 

S(w) = Fourier transform of the quadrature function S(t), 

fjs( W ) ] r » = o 

and S(w) = JO >> if -j w = 

[-jS(«)J U>o 

A 

Let s(t) = s(t) + j s(t); 

A A 

then s(t) ** S(w) + j S(w) = S(w), 



A fS(w) + j S(w) ] f w < 

S(w) - ^ S(w) ^ if J w = 

I "JS(w)J U> 

1S(w) [ 
US(»)J 



A [0 f u < 

S(w) -JS(w) [• if -j w = 

u> > 



40 



APPENDIX C.-SYMBOLS 

b(t) time-domain equivalent for boxcar function 

B(a>) frequency-domain equivalent for boxcar function 

e(t) envelope (instantaneous amplitude) 

f frequency 

f maximum frequency within wave train 

f,(t) instantaneous frequency 

f nq Nyquist frequency (folding or aliasing frequency) 

f(t) forcing function (signal) 

F(w) transform equivalent of forcing function 

F*(w) conjugate transformation of forcing function 

i(t) impulse response function (operator, filter) 

I(w) transform equivalent of impulse response function 
(transfer function) 

Im(c<;) imaginary part of Fourier spectrum 

\u>\j) cross-phase spectrum 

N number 

NP2 number of discrete points associated with time 
sequence 

o(t) transient response function 

O(w) transform equivalent of transient response 
function 

R(u>) real part of Fourier spectrum 

S(f) frequency-domain representation of analog signal 



USED 

Sgn(w) 

s(t) 

A 
s(t) 

S(«) 



IN THIS REPORT 

Signum function 

time series function (signal) 

Hilbert transform of s(t) 

frequency-domain equivalent of s(t), Fourier 
transform of s(t) 



A 
S(«) 


Fourier transform of s(t) 


|S(«)| 


amplitude spectrum of signal s(t) 


t 


time 


T 


period of window 


W c (t) 


cosine weighting function 


w h (t) 


Hanning weighting function 


WJt) 


half-sine weighting function 


t 


phase 


M) 


autocorrelation function 


*«(») 


frequency-domain representation of <^ ff (t) 


h(i) 


correlation function 


*■(») 


cross-power spectrum 


*i(0 


instantaneous phase 


4>u> 


phase spectrum (or angle) 


w 


angular frequency 



INT.BU.OF MINES,PGH.,PA 29088 



"0 

m 

z 

£ Q 



3 

0) 



ro 

2.2 

m 



03 cz 
Sco 



T1 

o 


> 


■o 


I - 


3 


00 


1 


c 


H 


CO 


m 

c 


z 
m 


m 


co 


I 


CO 



3 

<£ co 
° 3 

b z 
■c* co 

00 

o 
o 



ii a) 

It 



3 



3 



> 

z 

m 
D 

c 
> 

o 

"0 

o 

3J 



3 

m 



O 

-< 
m 

3J 



271 ^0 



^ 




o V 




.^Va\ ^ ^ .vote*. *x 



w 
**** 






Ik.X. <?.<$&>» .^.-'J^i'S. #+s2&.% >>..j&>X t o**i-.. 















.** /jgfe\ v>* .^va\ ^ .** /^fei # . %> J /aVa\ ^ ^ y^K*. ^. ** /; 



A 






A- • "« 






£ °* /^4r*\ c *.!^fc.% /^i^ 






^o 1 







V **TTV 



S 5 ^ 



%\ ^ im>y* ^s-^Sr <r time's V /^s%" c°*V^X 

^..Z "-o;-^S\o^ '° X**?P\/ °v^R % '.^ '°^v^? ; ./ ^o/^S'^o 










Av .HI. "*i 




^ ^ -.'« 









» "7 







->' 


















v 5 ^ 















^* -N v 









% 



•->•>*-■ 









V s^A>* ^ 








+**. 
■<?, 



**<? 



«V * • . o A 
















